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I. INTRODUCTION 



Determining the structure of the nucleon in terms of quarks and gluons is central to our goal of understanding 
baryonic matter at the level of its smallest constituents. While the theory describing the strong interactions of quarks 
and gluons, Quantum Chromodynamics, was identified thirty-five years ago, its predictions at low energies have been 
notoriously hard to derive ab initio. The modern approach to calculate the properties of hadrons is based on the 
Euclidean path integral representation of QCD discretized on a space-time lattice, i.e. lattice QCD. Importance 
sampling methods, implemented on massively parallel computers, make it possible to extract, in particular, many 
properties of the nucleon. 

In recent years, advances both in algorithms and in computer technology made a series of remarkable calculations 
possible that had a large impact on our understanding of nucleon structure. Among the quantities calculated we 
would like to mention the quark contribution to the nucleon spin ^ [2], the nucleon transverse structure [3], and 
the nucleon axial charge |4j |5] . Recently, the nucleon electromagnetic and axial form factors have received special 
attention in Refs. [B] and [7] using dynamical Wilson and asqtad fermions. Another important milestone is the advent 
of full domain- wall calculations, see Refs. [HUH], an d of dynamical twisted-mass fermions [TU]. Disconnected diagrams 
play a key role in an ongoing study of the strange quark content of the nucleon |11| . For reviews and progress reports 
on the current state of the field, see Refs. [T2H15] . 

Over the past years several of us have reported on hadron structure measurements using mixed action calculations 
with 2+1 flavors of dynamical asqtad sea quarks [T5], [T7] — corresponding to degenerate u,d quarks + the strange 
quark — and domain wall valence quarks |18H22j . A significant milestone was reached in Ref . |20j , which summarized 
our findings for higher moments of generalized form factors. The current paper represents a major update of that 
work: it includes the observables presented previously with higher statistics, as well as an additional, lower pion mass 
calculation. Beyond that, it covers form factors and chiral extrapolations of the forward moments that were not shown 
previously. The propagators and technology underlying these calculations have not only successfully been applied to 
nucleon structure, but have also turned out to be enormously valuable to other studies, see e.g. [2"5H2"5] and references 
therein. Mixed action calculations have also been studied in the framework of effective field theory, see 26 28J . It is 
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the purpose of the current paper to report our final results including all improvements we have made over the years 
on the method and technology of our computations. 

The layout of this paper is as follows. We give an overview of our notation and conventions and present the 
observables that we study in Sec. [TTJ The discussion of our technology and of the improvements we have made over 
our previous calculations together with consistency checks takes place in Sec. |III| Our results are presented in Sec. |IV| 



which is divided into the following subsections: Section IV A is dedicated to the axial charge qa- The discussion 



of the electromagnetic form factors and the axial form factors takes place in Sec. |IV 5] and [iVCl respectively. The 
generalized form factors and their extrapolation using different schemes of chiral effective theory are presented in 
Sec. HVDl 

The generalized form factors of the energy-momentum tensor provide vital information as to how spin is apportioned 
within the nucleon — a long-standing puzzle of hadron physics. In particular, they enable a first-principles calculation 
in lattice QCD of J q , the total angular momentum carried by the quarks [TJ [21 [Ml HS1 [50] , and hence have a crucial 



role in resolving this puzzle. Section IV E is dedicated to this topic. Our summary and outlook for future work are 
given in Sec. |Vj Since the summary includes cross-references to the most significant tables and figures, the reader 
mainly interested in new results might find it useful to use that section as a guide to the highlights of our calculations. 

II. OVERVIEW OF PHYSICAL OBSERVABLES 

The observables we are reporting on are defined via matrix elements of bi-local light cone quark operators in nucleon 
states. They can be systematically discussed in the framework of so-called generalized parton distributions [5TM34) . 
The relevant bi-local operator is given by 

O q , T {x)=J ^q^r^Ve-^^^^Ui^), (1) 

with x being the momentum fraction, n a light-cone vector, and T representing any gamma matrix from the basis 7^, 
7m 75, and a M , j = 1,2. 

The quark field q in Eq. ([I]) can carry any of the up, down or strange flavors, however here we restrict ourselves to 
q = u,d. Since in our lattice calculations the up- and down-quarks have degenerate masses, isospin symmetry is built 
in by construction. Note that our matrix elements always refer to the case of a proton, i.e. the nucleon sources used 



in the construction of the matrix elements in Eqs. ( 12 1 and ( 13 ) below always contain two u- and one e?-quark. When 
comparing lattice matrix elements to experiment, we need to choose appropriate flavors or flavor combinations. The 
isovector combination O q = O u — Od, where we subtract the down-quark contribution from that of the up-quarks, 
can be constructed from the difference between proton and neutron observables obtained in experiment. The isoscalar 
combination O u +d, on the other hand, corresponds to the sum of proton and neutron observables. To compute proton 
observables as measured by probing a proton with a photon, one needs to take into account the charge weighting 
factors of the quark-photon vertex and thus consider the combination q — 2/3u — l/3c?. In all but the isovector case, 
the matrix elements will in principle receive contributions from both the connected and the disconnected diagrams. 
In this work we have neglected the disconnected contributions since they are very costly to obtain. For recent studies 
of disconnected contributions in the framework of calculations of form factors and moments of parton distribution 
functions (PDFs) we refer to, e.g. Refs. [551 1551 . 

The matrix elements of the operator Eq. (fTl) between nucleon states with momentum p and polarization A can 
be parameterized generically by exploiting their Lorentz tensor structure. This has been discussed in detail in the 
literature [53] and we merely present the results here. 

For the case = 7^, the nucleon matrix element adopts the form 



(p',\'\O qn (x)\p,\) = ((fi))H(x,t,t) + ^K a »4^ 'i ■ ^ 



where we have introduced the notation ({X}} = u(p' , X')Xu(p, A) and the parameters A = p' —p, t = —Q 2 = —(p' —p) 2 



and £ = — n ■ A/2. In the framework of form factors in Sees. IV B and IV C we will denote the squared momentum 
transfer by Q 2 since this is a common and wide-spread convention. The unpolarized generalized parton distributions 
(GPDs) H(x,^,t) and E{x^,t) are Lorentz scalars and thus frame-independent functions parameterizing the matrix 
element. We point out that the matrix element also depends implicitly on a renormalization scale, /j 2 , and scheme. 

In the case T ^ = 7^75 we obtain a Lorentz-covariant parameterization in terms of the polarized GPDs H (cc, £, t) 
and E(x, £, t): 

77 • A 

(p',X'\O qni5 (x)\p,X) = ((fi l5 ))H(x,^t) + ^—(( l5 ))E(x,^t). (3) 

2m 
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The case r MJ = a^j will be discussed in a separate publication [37] . The kinematic parameter x is the average 
longitudinal momentum fraction of the struck quark and £ and t are the longitudinal and the total squared momentum 
transfer to the nucleon, respectively. The GPDs are defined over the full intervals x — — 1 • • • + 1, and £ = — 1 • • • + 1. 
Depending on whether |jc| > |£| or vice versa they have the interpretation of amplitudes for the emission and absorption 
of a quark or for the emission of a quark-antiquark pair, respectively. 

An attractive feature of the generalized parton distributions is that they occur in a range of different processes, 
e.g. deeply virtual Compton scattering, wide-angle Compton scattering and exclusive meson production, in addition 
to the classic processes that probe the forward parton distributions and form factors. The challenge of GPDs lies in 
their more complex structure — each generalized parton distribution is a function of three parameters rather than just 
one, and the different experimental processes provide different constraints on their form. Typically only convolutions 
of these functions in the x variable are experimentally accessible. 

Since lattice calculations deal with operators and matrix elements in Euclidean space, a direct computation of 
non-local light-cone elements is not possible. To facilitate the lattice calculations, one takes x n_1 -moments of Eqs. ^ 
and ([3]) , yielding a tower of local operators whose matrix elements can be related to the corresponding moments of 
H , E, H and E. In this study, we will compute matrix elements of the following local generalized currents, 

{Mi...M n } = g(o)r^ 1 iB M2 ---izWg(0), (4) 

where, again, r M can refer to either 7^ or 7^75. Curly braces around indices represent a symmetrization and the 

<-> — * <- 

subtraction of traces of the indices and the derivative is defined via D— 1/2(Z? — £))■ 
Taking the moments w.r.t. x of the GPDs we define 





J dxx n ~ 




E n (£,t) = 


J dxx n ~ 




EE 


J dxx n ~ 




= 


J dxx n ~ 


^E(x,^t). 



(5) 

The non- forward nucleon matrix elements of the local operators, Eq. Q, can in turn be parametrized according to 
their Lorentz structure in terms of generalized form factors (GFFs) A nm (t), A nm (t), B nm (t), B nm (t), and C nm {t), 

(p',X'\0^\p,X) = ((^))A w (t) + ^-((a^ a ))^ a B w (t), 

(p',X\0^^\p,X) = p^((r 2} ))A 20 (t) + ^-p^((a^ a ))A a B 20 (t) + -A^A^C 20 (t), 

2m 771 

(p',X'\0^^\p,X) = p { ^p^((Y l3} ))A 3Q (t) + ^-p^p^((a^ a ))A a B 30 (t) 

2m 

+A^A^((r 3} ))A 32 (t) + ^-A^A^((a^ a ))A a B 32 (t) , (6) 

2m 

for the vector operators and 

(p',X'\0%\p,X) = ((^ 5 ))A w (t) + ^-A^((To))B 10 (t) , 

(p',X'\0\^\p,X) = p{^((7^>7 5 ))i 20 (t) + J-A{"^><(76))B 2 o(t) , 
,a 2m 

(p',X'\0^^\p,X) = p { ^p^((r 3} l 5 ))A 3 o(t) + ^A^p^p^(( l5 ))B 30 (t) 

+A^A^((7^> 75 ))i 32 (t) + ^-A^A^A^(( l5 ))B 32 (t) , (7) 

2m 

for the axial vector operators. Here we have defined the average nucleon momentum p = (p' +p)/2. By comparing 
this with the x n_1 -moments of Eqs. |2]), and using Eq. (§, one finds that the £ dependence of the moments of 
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the GPDs is merely polynomial, 



H n, 




= A 10 {t), 




H n- 




= A 20 (t) 4 


(26) 2 C 2 o(t), 


H n- 




= A 30 (t) + 


(20 2 A 32 (t), 


E n, 




= B 10 (t), 




E n- 


- 2 (tt) 


= B 2Q (t)- 


(2£,) 2 C m (t), 


E n, 




= B3o(*)4 


(20 2 S 32 (<), 


H n - 


=1 (t,t) 


= A w {t), 




H n= 


= 2 (U) 


= A 20 (t), 




H n - 


- 3 (U) 


= A 30 {t) + 


m) 2 A 32 (t), 


E n, 


=1 (U) 


= B w [t), 




E n, 




= B 20 (t), 




E n, 


= 3 (U) 


= S 30 (i)4 


m 2 B 32 (t), 



(8) 

and 



(9) 

In the forward limit of Eqs. ^ and ^ with p — p', we obtain the well-known parton distribution functions, 

q(x) = H(x,£ = 0,t = 0), 

Aq(x) = H(x,£ = 0,t = 0). (10) 

Note that in the case £ = the GPDs — and also the GFFs, including the form factors for n = 1 — admit a 
probability interpretation [38] and that this property holds even in the case t ^ 0. Taking together Eqs. Q, ^ 



and ( 10 1 and setting t — will similarly yield the moments 

(x"- 1 ), = H n (0,0) = A n0 (0) , 
(a*- 1 ) A, = H n (0,0) = A nO (0). (11) 
The matrix elements are obtained on the lattice from the two point functions 

C^(T,p) = ^V 1 ** Tr (r poI (n(f, T)n@, 0))) , (12) 

and the three-point functions 

C 3 ^(T, T ,p,p") = £ e-V'-t+W-^ Tr (r pol (n(x, T )O(y, T)n(0, 0))) . (13) 

S,y 

We have introduced the lattice proton operators, n(x,T) and n(x,T). In order to maximize overlap with the ground 
state, we use the smeared sources defined in Ref. [39] . This overlap can be parameterized by a function Z(p) 
according to (Q\n(x)\p, A) = \J Z(p)u(p, \)e~ lp x . We also use the projection operator, r pol = |(1 + 74)5(1 — 17375)- 
O = q(0)Jq(Q) denotes the operator with all appropriate indices in which we are interested. 

Applying the transfer matrix formalism yields the following expression for the behavior of the two- and three-point 
functions 

Z(v)er ET 

C 2pt (T,p) = yFJ Tr(T pal (tf + m N )) + excited states. 

2E 

JZ(v)Z(v'V- e ' {Tq - t) - et 
C 3 *(T,T ,p,f) = 2E'2E Tr(r pol (i^' + TO JV )^(i^ + m Ar ))+ excited states. (14) 

In order to cancel the exponential factors and wave-function normalizations, we construct the ratio 



Ro{TTq) = C 3 o pt (T,T ,p,p>) j C^jTo - T,p)C^(T,p>) 



v /0*(T,£)O*(T > £') V ° 2pt ( T ° - TJ')C^{T,p) 

T »T», y u(p,\)r pol u(p' ,_X) (p 1 X'\Q\p X) , (15) 

y/2E(E + m N )2E'(E' + m N ) * 
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which becomes proportional to the desired matrix elements for sufficiently large source-sink separations, Tq, and with 
the operator insertion sufficiently far from both source and sink, T> 1 and To — T 3> 1. Typically, we find a plateau 
region [T min , T max ] over which we average the resulting value of the operator. 

Since we operate with a finite lattice of extent aL, the momentum values which we can choose are discrete and are 
given by p = 2n j (aL)n p with ri p a vector whose components are integers ranging from —L/2 to L/2. For the nucleon 
sink we choose the two values p' = 2ir/(aL)0 and p' = 2tt / (aL)(— 1, 0, 0) and for the source we choose p such that the 
absolute value of the integer momentum vector, \n p \ is smaller than five. This defines the set of t values accessible in 
our calculation. Note that the (generalized) form factors will receive contributions from several different momentum 
and index combinations at any fixed value of t. By constructing an overdetermined system of equations from all those 
combinations we make optimal use of the available data. This procedure has been discussed in detail in Ref. [2]. 

The energy of a state at momentum p is related to its mass through the dispersion relation. In our analysis we 
use the continuum dispersion relation. We have verified that the resulting energy agrees, for the spatial momenta 
employed in our calculation, with the energy of a nucleon at non-zero momentum actually calculated on the lattice. 

In the case of the electromagnetic current, ip^/jiifi, the generalized form factors correspond to the electromagnetic form 
factors of the nucleon. For the axial current, ip^^^ip, the form factors correspond to the axial and the pseudoscalar 
form factors. These will be covered in detail in Secs. lIVBl and lTYCI 



III. NEW LATTICE CALCULATIONS 



We now present the methods and technologies we have used for our calculations. As discussed previously in Ref. |20) . 
we continue to employ the asqtad action for the sea quarks and the domain wall (DWF) action for the valence quarks. 
In addition, we also add one lighter mass to our data set. The data sets used in this paper are summarized in Tab. [I] 
The columns show the bare asqtad quark mass for Nf = 2 + 1 dynamical fermions, the corresponding bare DWF 
mass, the volume in lattice units, the number of gauge field configurations used and the number of measurements 
included in the analysis. The bare valence DWF masses have been tuned such that the physical pion masses agree 
with those obtained from the purely asqtad calculation. The size of the 5th dimension has been set to L 5 = 16. The 
choice and tuning of these parameters has been discussed in detail in Ref. [2U] . 



m a,qtad f JVy =2+1 




Volume 


# confs 


# meas 


0.007/0.050 


0.0081 


20 3 


x 64 


463 


3704 


0.010/0.050 


0.0138 


28 3 


x 64 


274 


2192 


0.010/0.050 


0.0138 


20 3 


x 64 


631 


5048 


0.020/0.050 


0.0313 


20 3 


x 64 


486 


3888 


0.030/0.050 


0.0478 


20 3 


x 64 


563 


4504 


0.040/0.050 


0.0644 


20 3 


x 32 


350 


350 


0.050 


0.0810 


20 3 


x 32 


425 


425 



TABLE I: Summary of our data sets. 



We work in a mass independent scheme. The lattice spacing is therefore independent of the bare quark mass, and 
its value has been determined to be a = 0.1241(25)fm, corresponding to a -1 — 1.591(32)GeV in Ref. [2T], taken from 
heavy quark spectroscopy [JD]- This yields a physical volume of V = (aL) 3 = (2.5 fm) 3 on the L 3 — 20 3 lattices 
and of V = (aL) 3 = (3.5 fm) 3 on the L 3 — 28 3 lattice. The physical values of the nucleon masses, pion masses, and 
pion decay constants are needed for our computation of hadron structure. These have been previously determined in 
Ref. [3T]. We have listed them in Tab. [TT| The columns show the bare asqtad sea quark mass, the lattice size, and the 
resulting lattice pion mass, pion decay constant and nucleon mass. Finally, they are converted to physical values in 
MeV 1 . 

In contrast with our previous publication |20j . in this work we did not use Dirichlet boundary conditions on the 
first and the middle time-slice, forming so-called chopped lattices. Instead, the technology that has been employed 
consists in computing multiple source/sink pairs on a single gauge field configuration; it is discussed in more detail 
in the next section III A We find this technique both more convenient to use and more powerful in making optimal 
use of the existing resources. The quality of the results is superior since it allows us to process eight source/sink 



1 Note that unlike in Ref. 1211 . in this paper we use the normalization of such that /„- = 92.4MeV 



G 



Light m: o 5 f ad 


Volume fl 


(am), 


(a/)* 


(am)jv 


[MeV] 


U [MeV] 


m N [MeV] 


0.007 


20 3 x 64 


0.1842(7) 


0.0657(3) 


0.696(7) 


292.99(111) 


104.49(45) 


1107.1(111) 


0.010 


28 3 x 64 


0.2238(5) 


0.0681(2) 


0.726(5) 


355.98(80) 


108.31(34) 


1154.8(80) 


0.010 


20 3 x 64 


0.2238(5) 


0.0681(2) 


0.726(5) 


355.98(80) 


108.31(34) 


1154.8(80) 


0.020 


20 3 x 64 


0.3113(4) 


0.0725(1) 


0.810(5) 


495.15(64) 


115.40(23) 


1288.4(80) 


0.030 


20 3 x 64 


0.3752(5) 


0.0761(2) 


0.878(5) 


596.79(80) 


121.02(34) 


1396.5(80) 


0.040 


20 3 x 32 


0.4325(12) 


0.0800(5) 


0.941(6) 


687.94(191) 


127.21(78) 


1496.8(95) 


0.050 


20 3 x 32 


0.4767(10) 


0.0822(4) 


0.991(5) 


758.24(159) 


130.70(67) 


1576.3(80) 



TABLE II: Hadron masses and decay constants in physical units. 



pairs instead of just two as before. We have recalculated our results on the lattices with asqtad sea quark masses of 
TO aaqtad _ 0.007-0.030. The higher masses are included in some plots to guide the eye, but have never been included 
in the chiral fits. Section [ill B| discusses possible systematic errors of our nucleon matrix elements. 

We also take all possible sources of correlations into account carefully by performing all fits using the error correlation 
matrix among all data points at fixed pion mass and the "super jackknife" technique for combining data from different 
pion masses in a single fit. These techniques are discussed in Sec. |III C| Section |III D| specifies our renormalization 
procedure for the lattice operators we use and lists all relevant renormalization constants. Finally, Sec. |III E| discusses 
potential influences of finite-volume effects on our lattice data. 

A. Improved statistics 

In order to improve the statistical quality of our data set, we employ a method we call "coherent sink technique". 
This method proves to be a substantial improvement over previous methods employed for the extraction of three-point 
functions from lattice data. Effectively, we obtain eight measurements of matrix elements per configuration. Hence, 
this allows us to make optimal use of the rather expensive configurations generated with dynamical quarks. 

On every other configuration we place sources at space-time positions (0, 0), (L/2, 16), (0, 32), and (L/2, 48), where 
L denotes a spatial vector with components (L,L,L). We then perform twelve inversions of the Dirac operator, 
corresponding to the four spin and three color indices of the quarks, and feed them into the construction of the 
forward propagators. Using these forward propagators, we create a momentum projected nucleon sink a temporal 
distance To away from the source, i.e. we end up with four sinks located at Euclidean times To, To + 16, T) + 32, and 
T + 48. 

If we did a conventional calculation, we would consider each source-sink pair completely separately, constructing a 
set of backward propagators from the sink and evaluating the three-point function. Thus, each measurement would 
require a separate set of backward propagator calculations. With our new approach, however, we calculate a single set 
of coherent backward propagators in the simultaneous presence of all four sinks! Combining these coherent backward 
propagators with the forward propagators yields the physical matrix elements plus terms that vanish due to gauge 
invariance when computing expectation values. 

Additionally, we also create a momentum projected antinucleon sink a temporal distance — T away from each source 
and perform an analogous calculation for coherent antinucleon propagators. It is then straightforward to relate the 
matrix elements of our twist-two quark operators in an antinucleon to the desired results in a nucleon. To summarize, 
given a set of forward propagators, we obtain eight measurements of Cq* at the cost of two rather than eight sets of 
inversions for the backward propagators. To minimize correlations, we alter the source locations on every other lattice 





coherent antiproton sinks 






coherent antiproton sinks 






























































































1 coherent proton sinks 






| coherent proton sinks 





FIG. 1: Layout of smeared nucleon sources (stars) and coherent, fixed momentum sinks (vertical lines) on our lattices: (left) 
even configurations; (right) odd configurations. 
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to be (L/2,0), (0, 16), (L/2,32), and (0,48). Thus, each source is shifted by a displacement of L/2. This layout of 
the sources is illustrated in Fig. [I] 

We find that the individual measurements only exhibit minimal correlations and thus this strategy provides a 
valuable increase of statistics, see Ref. [U] for a discussion of autocorrelations. We consistently apply binning with a 
bin-size above ten which eliminates residual autocorrelations from our data. Below we study whether this technique 
introduces any systematic errors into our calculation. 



B. Tests for systematic errors 



We have performed tests concerning three possible sources of systematic errors: The coherent sink scheme described 
in the previous section, the choice of boundary conditions in the temporal direction, and the choice of the source-sink 
separation, To, which is an important input parameter entering any three-point function calculation 2 . Following 
Eq. ( 14 1 it is advisable to pick To as large as possible. However, if Tq becomes too large the two-point function in 



the ratio Eq. ( 15 ) introduces an increasing noise that will eventually wipe out the signal. Thus, To should be chosen 
such that we are still able to project out the ground-state in a suitably chosen plateau-region, but not so large that 
the signal to noise ratio becomes too bad. 

From previous calculations, cf. e.g. [50], we know that a separation of To ~ 1.2 fm is a reasonable choice. Our tests 
of this assumption, together with the other characteristics of our calculation mentioned above, are summarized in 
Fig. [2j The labels for the calculations denoted on the abscissa are explained in Tab. |III| The gauge fields are a sample 
of 448 configurations at the working point ™^ tad = 0.010/0.050 on the 20 3 lattice, cf. Tab.Jl] The technique employed 



1.5 
1.4 
1.3 
1.2 
1.1 
1 



BCD 
Calculation 



25 r 



22.5- 



J 20 



17.5- 



15 L 



BCD 
Calculation 



FIG. 2: Comparison of coherent sink technique and tw o source-sink separations with our previous calculation. The obser vables 
are the axial charge, gu in the left panel, see Sec. IV A and the isovector charge radius, a~ 2 {r\) in the right panel, see Sec. IV B 



Label 


# meas 


To 


Technology 


A 


448 ■ 8 


9 


Independent backward propagators, nucleon/antinucleon 


B 


448 • 8 


9 


Coherent backward propagators, nucleon/antinucleon 


C 


448 ■ 8 


10 


Coherent backward propagators, nucleon/antinucleon 


D 


448 


10 


Single source/conf., nucleon only, unchopped L t = 64 lattice 


E 


448 


10 


Single source/conf., nucleon only, chopped L t = 32 lattice 



TABLE III: Calculation techniques used for the comparison of the coherent sink techniques and two source-sink separations 
with our previous calculation. 



in our previous publication, Ref. [20], corresponds to the label "E", cf. Tab. Ill To address whether the chopping 
prescription — to cut the lattice into two halves, impose Dirichlet boundary conditions, and compute observables on 



2 The issue of finite size effects is addressed separately in Sec. 
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both halves separately — employed in that publication introduces a systematic error, we have repeated the calculation 
on the unchopped lattice. This case is denoted by label "D". It is evident that the results are essentially identical 
both for the axial charge and the mean squared radius. We thus conclude that the chopping prescription did not 
introduce any systematic error. 

Next, using the same sample of configurations we performed eight independent calculations of propagators on the 



lattice with locations identical to those chosen for the coherent propagators, cf. Sec. Ill A This case is denoted by "A" . 
It is evident that the statistical improvement is remarkable, indicating that the results are sufficiently decorrelated to 
warrant the extra effort. 

The true power of the coherent sink technique is demonstrated by case "B" . The computational requirements for the 
backward propagators are just 1 / 4th of those of case "A" , yet the result is almost identical and completely consistent. 
This gives us further confidence that our technology is indeed correct. 

Finally, we also compare the two source-sink separations, To = 9 and To = 10, the latter denoted by case "C" . There 
is an increase of about 50% in the error bar between case "B" and case "C" , but the results are fully compatible within 
error bars. This indicates that a separation of To = 9 is already sufficient to extract the ground state and higher-state 
contaminations are negligibly small. We thus proceed to use the method denoted by "B" for the remainder of this 
publication. 

To illustrate the result from different source-sink separations and to make sure that excited state contributions 
are small, we also show a plateau plot for one of our observables. Figure [3] shows the isovector axial charge, gA, as 
a function of the location of the operator insertion prior to averaging over the plateau region. The plateau is flat 
up to 0{a 2 ) cutoff effects for any To since the operator corresponds to a conserved charge. If we had to suspect 
contamination from excited state contributions, the overall level of the plateau would be systematically different from 
the corresponding ones using the To = 9 separation. We do not observe such a systematic effect and thus conclude 
that going to the smaller separation does not introduce systematic uncertainties from excited states. 
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FIG. 3: Comparison of two separations on two volumes for the isovector axial charge, gA- This plot shows the resulting plateau 
plot, i.e. <m(T) as a function of the location of the operator insertion at fixed source-sink separation. 



In summary, the new technology does indeed constitute a major improvement over the one previously employed in 
terms of statistical accuracy, but does not lead to a detectable increase in systematic error. 



C. Super jackknife analysis and error correlation 

When calculating physical observables and their uncertainties from different lattice data sets we have to take into 
account the correlations of different data points with each other. It turns out that there are three possible cases 
we need to consider: (a) Observables calculated from data points computed on a single ensemble of lattice data, 
(b) Observables calculated from data points computed on several statistically independent ensembles of lattice data, 
where each ensemble contributes only a single data point, (c) Observables calculated from data points from different 
ensembles, where each ensemble contributes more than one data point. An example for case (a) are dipole fits to form 
factors, where each fit is done to a set of different Q 2 values at a fixed pion mass. Case (b) occurs for chiral fits to 
hadron masses [21 or moments of forward parton distributions [39j . Case (c) is the most complicated and shows up 
in simultaneous fits to (generalized) form factors as functions of both Q 2 (or of t in case of the GFFs) and m n . 

In case (a) we are dealing with data that may have correlations among data points. Case (b) has no correlation 
among different data points since the underlying gauge field configurations are entirely independent. Case (c) has data 
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points that are correlated (those points obtained on the same sample of gauge field configurations) and others that 
are not correlated (those obtained on different samples of configurations). Although in case (c) the error correlation 
matrix will be strictly block-diagonal, a straightforward numerical estimation may not take this property into account. 
In particular, the off-diagonal entries in the error correlation matrix may have large uncertainties themselves and thus 
introduce numerical instabilities. 

In this paper, we have decided to consistently adopt the jackknife method to compute uncertainties and to use the 
error correlation matrix in the function x 2 where it does not introduce numerical instabilities. The standard jackknife 
error prescription is well-known and discussed widely in the literature, see e.g. |42l 143] . We briefly summarize the 
method for an observable with exact statistical mean A as follows: given a series of measurements {a^}, i = 1 . . . N, 
with N being the total number of measurements, dj being the i th measurement we define the i th jackknife average via 

1 

(16) 



N-l ^ 



Thus, we obtain a new set of N jackknife averages, {a.;}. The natural estimator for A reads a — X)i=i ®i m terms of 
the {a-i}. More generally, for any function f(A) of the observable, the set of N jackknife blocks provides an estimate 
of the mean F, 



i=l 



as well as an estimator for its uncertainty, 

1 ]T(/(a,)-/(«)) 2 . (18) 



N-l N 



i=i 



Note that this method also generalizes to the case of a function of several observables, labeled by a Greek index, with 
expectation values A a , a — 1 . . . n. Again, this topic has been widely discussed in the literature, including in our 
previous studies, see e.g. (25]. As an important example, we describe the case of a fit to several observables. The 
function that we minimize in a nonlinear fit to n different lattice measurements is given by 

n n 

^ = E !>«({*}) - - *p,i)C~l , (19) 

o=l /3=1 

where y a ({x}) is a model function with parameters {x} and a possible dependence on the index, a. The covariance 
matrix C a p is defined by 

JV-1 N 

C a l3= E(° a ~ a <xj)( a P ~ £.*) ' ( 20 ) 

i=l 

For each configuration with number i we substitute the jackknife averages a a j and minimize xf w.r.t. the parameters 



{x}. The parameter values which minimize Eq. (19) are thus implicit functions of the original data set. Finally, 
we use Eqs. (17) and (18) on the resulting set of parameters to obtain estimates of their central values and of their 
statistical uncertainties. 

As a special case, we consider the case of a dipole fit at fixed pion mass, where each a a ^ denotes the form factor 
measurement on the i th jackknife block at momentum transfer Q^, and the function y a (Md, Aq) = Aq/(1 + Qa/M^) 2 



in Eq. ( 19 1 is the dipole model function, cf. Eq. ([37]) below, the parameters being the dipole mass, Md, and the overall 
normalization, Aq. The best fit parameters are obtained from a minimization of the associated function % 2 for all 



jackknife blocks in the sample and we get the results by using Eqs. (17) and (18), where the function f({A a }) denotes 
either Md or Aq. 

As has been mentioned previously, it may happen that the off-diagonal elements of C a p are only poorly determined 
and the resulting inverse has zero or negative eigenvalues. In that case, the function x 2 is not positive-definite and 
the problem is ill-defined. If this happens, we have to resort to simply using the diagonal matrix elements of C a p and 
the method reduces to the regular, uncorrelated fit. We point out that even in this case the statistical error associated 



with this procedure is correctly estimated by the jackknife method, Eq. (18). 

We also want to point out that each measurement can be taken from a single gauge field configuration or obtained 
as an average of a small set of B subsequent configurations. In the latter case, the jackknife procedure corresponds to 



10 




FIG. 4: Comparison of bootstrap (left panels) and super jackknife (right panels) resampling plans for chiral fits to the nucleon 
axial charge, <?a- The upper two plots show fits to the lowest three masses, the lower two show fits to the lowest four masses. 



what is called "jackknife with block-size B" in Ref. [35]. To summarize, the jackknife method allows us to compute 
functions of averages obtained on a given ensemble. The estimation of the uncertainty in Eq. ( 18 ) correctly takes into 
account correlations among the input data. 

When we are dealing with a function of several observables computed on different ensembles, the method needs to 
be generalized further. In order to take the particular features of such a "super-sample" containing both correlated 
and uncorrelated data into account, we use an extension of the jackknife method, called "super jackknife" [44l 145] . 
The idea is to define generalized (super) jackknife blocks, with which averages and errors can still be estimated using 
Eqs. (17) and (18). To illustrate this procedure we start with M distinct, uncorrelated ensembles with Nk samples 

available in the k th ensemble. On these ensembles we again have a set of measurements, {a^)}, which denotes the 



measurement on the i 



th 



sample of the k th ensemble. The averages affl = l/A^X^t a< ai an d the jackknife blocks 



,JV* (k) 



a a \ are introduced as in Eq. (16) above. The total number of super-jackknife blocks is defined by TV 



With appropriate super-jackknife blocks {5,^} we can generalize Eq. (19) to 



E E E Oi fe) (M) - *S) (»? } ({*» - ag>) {c^)-l 



(21) 



k=l a=l /3=1 



where we introduce the notation for the number of observables in the k th ensemble. The index i now denotes the 
number of the super-jackknife block, i = 1 ... N. The N super-jackknife blocks a^\, i.e. the N sets of arguments at 
which the function is to be evaluated are constructed as follows. The first TVi blocks consist of 



~{k) 



:* = 1 

ai fc) : k ^ 1 



with i = 1 . . . Ni, k = 1 . . . M, a = 1 . . . n 



(fc) 



(22) 
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The following N2 blocks consist of 

n [k) - =) ' with i = l...N 2 ,k = 1.. .M,a = l...n {k) , (23) 

aL fe) :k^2 

and so on. This generalization takes the correlations within each ensemble correctly into account, while at the same 
time implicitly sets correlations among different ensembles to zero. It is evident that it reduces to the regular jackknifc 



method, Eq. (16), in the case of a single ensemble. 

A typical case that the super-jackknife method can be applied to is the fit of a form factor which has been expanded 
simultaneously in and Q 2 to all available lattice data, see e.g. the SSE formula in Eq. (43). The function 



Va\cA, B[ (X)) — u(ca, BJ* (A);mi , (Qa ) 2 ) in Eq. (21) is the form factor model function and Eqs. (22) and (23) 
collect all form factor data on the ensembles with different pion masses, and at different (Q^) 2 values. Note 

that the values of Q 2 on each ensemble are different in general since both the lattice volumes and the nucleon masses 
are different. After N minimizations to the super samples, we obtain all parameters — in our example the low-energy 



constants ca and B£ (A) — by virtue of Eqs. (17) and (18) from the resulting set of best fit parameters 



In order to ascertain that results obtained with the super jackknife method are compatible with those obtained 
with competitor schemes, like the bootstrap resampling plan, which has, e.g., been employed in Ref. |21] . we have 
made a detailed comparison between the two methods based on fits to our data for the nucleon axial charge, qa- 
Figure [4] shows a comparison of two fits to qa-, one based on the three smallest pion masses and one based on the four 
smallest pion masses, corresponding respectively to the upper panels and the lower panels of the figure. The two left 
plots show error bands determined from the bootstrap resampling method, and the two right plots show results from 
the super jackknife prescription. The technical details and results will be discussed later in Sec. IV A| we use this fit 



merely as a test to verify that the bootstrap and the super jackknifc methods give indistinguishable error bands in 
the two situations. We point out that the case with three masses is a "bad" fit with large uncertainties and the case 
with four masses has substantially smaller uncertainties and provides a much better fit. We conclude that the two 
resampling schemes give essentially identical error estimates and are thus equally applicable, both in fits with large 
and in fits with small uncertainties. 

The error correlation matrix also allows us to study a phenomenon about lattice data that looks puzzling when 
one assumes that data points at distinct momentum transfers Q 2 are independent. In this paper we will plot a 



couple of those cases, see e.g. Fig. 25 A couple of data points are systematically higher, they are off by more than 
one standard deviation from the central fit. If these points were independent measurements, it would be a highly 
significant deviation. Thus, the question is whether these six points are highly correlated so that only one or two 
degrees of freedom have fluctuated randomly. Using the technology of error correlations presented in this section we 
have addressed this phenomenon in detail in Ref. [41] . We have found that the correlation of such data points is very 
high, typically (80 — 90)%. On the other hand, the resulting x 2 /dof is still around 1 which proves that such "outliers" 
are not statistically significant, although a naive visual analysis would lead to the erroneous conclusion that the data 
is incompatible with the fit function. In this way the treatment of error correlations as implemented in this work is 
absolutely necessary to derive conclusive statements about the statistical quality of fits to lattice data points. 

D. Renormalization of lattice operators 

The matrix elements we compute numerically are in the lattice regularization at a scale of the lattice cutoff, 
/i 2 = a~ 2 . In order to compare them to experiment, we need to convert them to a commonly used scheme. When the 
renormalization is multiplicative, the conversion is done via 

{O cont ) = Z {O utt ) , (24) 

(O cont ) being the matrix element in the MS-scheme and (0 latt ) the bare lattice operator. The factor Zq is the 
renormalization constant which depends on the details of the lattice action and the operator, but not the external 
states of the matrix element. 

Since on the lattice all operators are representations of the finite hypercubic group H(4), there will necessarily be 
fewer operators than in the continuum group, 0(4). This means that in general, continuum operators correspond 



to linear combinations of a finite set of lattice operators and the r.h.s. of Eq. (24) is replaced by a sum containing 
operators of different dimensions. In such a situation, computing renormalization coefficients will be impractical since 
subtracting power-law divergences non-perturbatively is required. It is possible, however, to compute the coefficients 
if operators of different dimensions do not mix — in the case of the operators appearing in Eq. Q, this requires 
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choosing distinct indices, which is only possible for operators with at most four indices. Hence, our lattice technology 
limits us to the computation of only the lowest moments of (generalized) parton distributions. 

We apply the perturbative renormalization constants computed in Ref. |46j . We have employed them previously 
in [5D]. All the perturbative coefficients relevant to this work are listed in Tab. IV together with the representation 
of H(A). In addition, we also obtain renormalization constants for the electromagnetic and axial currents non- 
perturbatively, see Tab.|V]in Sec. |IVB| and Tab. |XXXHl| in Sec. |IVC| In particular, the wave function renormalization 
encoded in the axial current renormalization factor Za implicitly enters most of the operators we study, and is not 
small; therefore, it is desirable to determine this one common factor non-perturbatively and employ the same recipe 
as in [20 , namely, 



— 7? ' "A, nonpcrt j 



Z A 



pert 



(25) 



with import = 0.964 for all but the vector and axial currents. The numerical results in this paper are all transformed 
to the MS-scheme at the scale of p? = 4GeV 2 . 



Operator 


H(i) 


/7Pcrt 


Q[ls]l{^D u} q 




0.962 






0.968 




t<*> 


0.980 


ghs]j{^D v D p} q 




0.982 



TABLE IV: Perturbative renormalization constants to convert bare matrix elements to the scale fi - 



in the MS-scheme. 



We note that an ongoing analysis of non-perturbative renormalization constants employing the Rome-Southampton- 
scheme [UJ turns out to be challenging due to discretization effects and the restricted range of applicability of the 
perturbative scale evolution. From our preliminary results [48] , we can however obtain a numerical estimate of the 
potential systematic uncertainty due to the renormalization of the one-derivative operators, which is ~ 7%. We 



include this systematic uncertainty in our final results in Table XLIV 



E. Finite-volume effects 



It is important to assess the potential influence of finite-volume effects on our observables. Our calculations are 
done in a fixed physical volume, hence we expect the finite-volume effects to increase as the pion mass decreases. 
Since at the working point m^f ad = 0.010/0.050 we have two volumes at our disposal, 20 3 and 28 3 , corresponding 
respectively to (2.5 fm) 3 and (3.5 fm) 3 , we have a way to estimate the size of these effects at a fixed pion mass of 
356MeV. 

We start with the simplest quantity, namely the mass of the nucleon. Figure [5] displays the mass plateau on the 
two volumes, together with the band from fitting the two-point correlators with the function [9 

C 2p t(i) = Zo mNt + Z ie - m ^ + (-lyz^e-z™* (26) 

in range 2 < (t/a) < 15. The same source is used on both volumes, and the effective mass is extremely similar. Fits 
to the mass plateau yield the result 

J7V=356MeV: |1 - m N (L = 2.bim)/m N {L = 3.5fm)| < 1.8% (95% conf. lev.), (27) 

As a consequence, in extrapolations carried out throughout this paper, we will use a common value of the nucleon mass 
for both volumes. Secondly, Fig. [3j used earlier to illustrate the source-sink separation dependence, also compares 
the plateau plot for gA on the 20 3 and the 28 3 lattices. The whole function qa (T), where T is the source-operator 
separation, is strikingly similar between the two volumes. We will return to this fact in Sec. IV A| dedicated to Qa- 



Figure shows the isovector Dirac radius, (r 2 ) extracted from a dipole fit to the form factor Fi(Q 2 ), as a function 
of the pion mass, m„. These observables will be discussed in detail in Sec. |IVB| here we only wish to exhibit the 
finite size effects. The upper cut-off in Q 2 for the dipole fit has been varied and the results for different cut-offs have 
been drawn with a slight displacement for clarity. At the pion mass of m w = 356MeV, there is a discrepancy outside 
the error bars between the two volumes. However, when reducing the data set to Q 2 values below 0.4GeV 2 , we find 
that the two Dirac radii are actually compatible. The discrepancy only becomes apparent at data points beyond 
Q 2 > 0.4GeV 2 . 
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FIG. 5: Mass plateau at am = 0.010 on the 20 and the 28 lattices. The band corresponds to Eq. ( 26 \ 
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FIG. 6: Isovector Dirac radius (rf } as a function of the pion mass for all ensembles. 



When studying the isovector form factor F^Q 2 ), we obtain the Pauli radii, (r|), and the anomalous magnetic 
moments, k v , shown in Fig. [7] Also these observables will be discussed in detail in Sec. IV B here we again only wish 
to exhibit the finite size effects. As in the previous plot, we have applied dipole fits with varying upper cut-off in Q 2 . 
Again, we find that there is a notable discrepancy for (r|) when data points at Q 2 > 0.4GeV 2 are included in the 
plot. For k v we find that the results on the two volumes deviate systematically, but are still compatible taking into 
account their statistical uncertainty. 

We can also scrutinize this behavior by looking at the form factors as a function of Q 2 directly and compare the 
location of the points on the two volumes. Figure [8] shows this comparison. The solid curves are dipole fits to all 
available data points. It is evident that their dipole masses are different, but the data points at small Q 2 are identical. 
The difference in curvature is only caused by points beyond that. 

We conclude that we observe finite-size effects in our lattice form-factor data. However, the finite-size effects are 
only significant for intermediate values of Q 2 . Based on the smallest values of Q 2 , the Dirac and Pauli radii may be 
identical on both volumes. We did not find a satisfactory explanation based on chiral expansions or models at this 
point and leave this matter for future investigations. 

As a final remark, in Sec. |IVD| we will compare the radii defined by the generalized form factors. There we will 
also find suggestive evidence that the mass radius is larger on the larger volume. Here our intention was to illustrate 
our observation of finite-volume effects with the data that is most accurate. 
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FIG. 7: Isovector Pauli radius, (rf), and anomalous magnetic moment, 
for all ensembles. 
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FIG. 8: Direct comparison of finite-volume effects for the isovector form factors, F^{Q 2 ) and F% (Q 2 ), at m n = 356MeV on the 
two volumes Q, = 28 3 x 64 and Q — 20 3 x 64. 



IV. RESULTS 



A. The axial charge 



In this section we present new data on the nucleon axial charge gA- Although axial form factors are discussed in more 
depth below in Sec. |IVC| the fundamental phenomenological importance of the axial charge warrants highlighting our 
new results for this observable already at this point. Aside from the changes in technology described in Sec. |III| our 
calculation follows closely the methods of Ref . [3] . In particular we use the local axial current for the calculation, and 
the five-dimensional axial current is used to determine the normalization of the local current, as described in detail 
at the beginning of Sec. |IVC| 

The new data are displayed in Fig. [9] The value of gA is remarkably independent of the pion mass, and lies at a 
value (8 — 10)% lower than the experimental value of 1.2695(29), while statistical errors are less than 2%. A naive 
extrapolation linear in of the < 500MeV data leads to gAim^) — 1.153(28). We will discuss below what 
difference more sophisticated chiral effective theory fits make. 

It is worth describing to what extent the situation has changed since the calculation [1] . In the latter, less accurate 
calculation, the lattice data also showed a very mild pion-mass dependence. Using a 3-parameter fit based on the 
leading one-loop pion mass dependence in the Small Scale Expansion (SSE) at finite volume leads to the value 
gA{rn n = 140MeV) = 1.226(84). The finite-volume effects predicted by the formula at the simulation points were 
found to be negligible compared to the statistical errors. The largest pion mass included in the fit was 760MeV, and 
the lightest 356MeV. 

Thanks to the new, higher statistics data, we control the finite-volume effects to a higher level of accuracy. Indeed, 
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FIG. 9: SSE fit to the axial charge. Left: three-parameter fit with < 500MeV. Right: two-parameter fit with < 360MeV. 



fitting the 20 3 and 28 3 g A plateau at am — 0.010, see Fig. [5J leads to the bound 

m ff = 356MeV : \g A (L = 3.5fm) - g A (L = 2.5fm)| < 0.045 (95% conf. lev.). 



(28) 



To further tighten this statement, we want to constrain the possibility that the plateau for g A could be affected 
by different excited state contributions on the two volumes. Indeed, even if the nucleon mass has a weak volume 
dependence for L > 2.5fm, see Eq. (27), the energy of the first excited state in that symmetry channel could a priori 



have a significant volume dependence: in large volume we expect it to be a nucleon and a pion with a non-vanishing 
relative momentum. However, comparing the local effective mass on the 20 3 and 28 3 lattices, Fig. [5J we see good 
agreement between them (the same nucleon interpolating operator is used on both volumes). We conclude that the 



contamination of the first excited state does not increase significantly with the volume. In particular, the bound (28) 



on the g A finite-size effects is robust. While a non-monotonic volume-dependence of g A that would make the difference 
in Eq. (28 1 accidentally small cannot be excluded, this bound strongly constrains how much of the discrepancy between 



the lattice data and the experimental value of g A can be attributed to finite-volume effects. 

As mentioned above, a naive extrapolation linear in m 2 of our lattice data leads to values of g^m^) about 10% 
lower than the phenomenological value. We now proceed with the fit ansatz provided by the small-scale expansion 
(SSE) framework 09], 
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(29) 



with gi the axial-delta-delta coupling, c A the axial-nucleon-delta coupling and A denoting the delta-nucleon mass 
splitting, in the chiral limit. Following |49j . we define the function 



R{m) 



A 

771 



A 2 



- 1. 



(30) 



When the A baryon is below threshold, as is the case in our lattice calculations, -y^A 2 
by ~\fn 



value of c A 



A 2 arccos(A/m7r) 
= 1.5 



m 2 log R(m„) is substituted 

A three-parameter SSE fit to our data at pion masses below 500MeV with a fixed 



a motivation for this choice is given in Ref. |S] — gives, see left panel of Fig. |9j 

1.22(17), .91=3.9(3.0). 



9° A 



(31) 



The result thus does increase as the pion mass is lowered, but only becomes consistent with the phenomenological 
value by virtue of its uncertainty also rising significantly. As an alternative, we perform a two-parameter fit to 
TOjr < 360MeV, where we fix the value of g\ to 2.5, close to the SU(4) spin-flavor quark symmetry prediction 9/5g A . 
This fit is illustrated on the right panel of Fig. [9j Here the result for g A is slightly lower than in the three-parameter 
fit, and the error bar barely extends to the phenomenological value. 
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As has already been mentioned above, our previous calculation of gA in Ref. [I] included data points at larger pion 
masses in the fit — as a result of which it has a smaller statistical uncertainty for the extrapolated value. Another 
calculation in Ref. i5j contains only pion masses larger than 500MeV. A discussion of the range of validity of one chiral 
expansion scheme in Ref. [50] concludes that lattice data below pion masses of 300MeV are necessary for a reliable 
prediction. Reference [ST] observes a bending down of the extrapolation due to the data point at the smallest available 
pion mass in that calculation, = 331MeV. This is above our smallest mass. However, our smallest pion mass data 
point has a larger error bar and is just consistent with the one from Ref. |51j : it may be the case that at this parameter 
the data is already affected by finite-volume effects - a possibility also mentioned in that paper. This interpretation is 
supported by the observation that the other data points at larger pion masses tend to be systematically higher than 
our data points. However, at our data point at m„ = 356MeV we do not find any evidence of finite-volume effects 
which indicates that at lighter pion mass these effects would have to set in rather quickly. 

While in the present work we find no significant evidence for a pion mass dependence of gA, our data is simultaneously 
compatible with the possibility that the functional form predicted by the small-scale expansion applies below — 
350MeV and with the phenomenological value of gA- 

B. Electromagnetic form factors 

The matrix element of the electromagnetic current between nucleon states can be parameterized in terms of two 
form factors. Common choices are the Dirac and Pauli form factors, Fx(Q 2 ) and F 2 (Q 2 ), and the electric and magnetic 
Sachs form factors, Ge{Q 2 ) and Gm{Q 2 )- The former directly correspond to the form factors ^4io(Q 2 ) and B 10 {Q 2 ) 
from Eq. (f6j). The latter are related by a simple linear transformation to the isovector Dirac and Pauli form factors, 
F\{Q 2 ) and F%(Q 2 ): 

G E {Q 2 ) = F v 1 {Q 2 )--^- 2 F^{Q 2 ) (32) 

G M {Q 2 ) = Fl{Q 2 )+F l 2 {Q 2 ). (33) 
We will also use the standard notation for the anomalous magnetic moment of the nucleon in units of e/2mN(Tn w ), 

Kv =F^(0). (34) 



When performing chiral fits we will work with 
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x - (35) 



m N (m^) 



which represents the isovector anomalous magnetic moment in units of the physical Bohr magneton, e/2rn^ lys . 

We use the ultra-local discretizations of the dimension three quark bilinear operators, i.e. their support is a single 
lattice site. Due to quantum effects the matrix elements of these lattice operators are not trivially renormalized, and we 
have to apply renormalization constants to them. Since the forward matrix element (p, X^j^iplp, A) counts the total 
number of quarks of type ip and this number is known by construction, we obtain Zy by dividing the unrenormalized 
isovector current in the forward case. We point out that in the forward case the disconnected contribution is exactly 
zero since the disconnected operator cannot change the total number of quarks of any type. Thus, the value for 
Zy obtained this way will be exact also if we consider disconnected contributions in future work. The resulting 
renormalization constants Zy for the vector current are listed in Tab. |Vj The renormalization constants of the axial 
current are discussed later in Sec. IIV CI 



m »sqtad 


Volume 


Z v 


0.007/0.050 


20 3 


x 64 


1.1159 


0.010/0.050 


28 3 


x 64 


1.1169 


0.010/0.050 


20 3 


x 64 


1.1206 


0.020/0.050 


20 3 


x 64 


1.1351 


0.030/0.050 


20 3 


x 64 


1.1464 



TABLE V: Renormalization constant of the vector currents. 



To study the charge distribution of the nucleon at large distances, it makes sense to consider the leading contribution 
of the form factors at small values of Q 2 [3J. The linear coefficient of the small-Q 2 expansion can serve as a measure 
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of the nucleon size and is known as the mean squared radius, 
structures one may consider: 



?), where i labels the different Lorentz and flavor 



Fi(0) ( l '//-' 



0(Q 4 ] 



(36) 



The radii, (rf), can also be extracted from experiment. For a recent review see Ref. 
forward for the proton isovector F^{Q 2 ) form factor, a determination from fits to the experiment 



Although this is straight- 
[SH turns out 



to be inconsistent with an analysis based on dispersion theory [55H57] . The latter radii are systematically larger 
than the former. To resolve this discrepancy, a dedicated experiment is currently being performed [58] . For the 
proton isovector F^iQ 2 ) a different discrepancy has been found in recent spin-transfer measurements |59H63j . The 
source of this mismatch is generally believed to be two-photon exchange processes [53] , which is challenging to verify. 
On the lattice, we can study these observables without any two-photon contamination and thus make a significant 
contribution towards resolving the discrepancy. 

This chapter discusses the res ults for the form factors of the electromagnetic current. First, w e study the isovector 



Dirac form factor F^(Q 2 ) in Sec. 



IVB1 



and the isovect or Paul i form fa ctor F: ? { Q 2 ) in Sec 



IVB2 



The scaling behavior 



of form factors at lar ger valu es of Q z is shown in Sec. IV B 3 Section IV B 4 discusses the Sachs parameterization of 



form factors. Section IV B 5 discusses the slope of the ratio Ff/Fi(Q 2 ) to learn about the flavor dependence of the 
form factors. The isoscalar form factors are shown in Sec. |IVB~6| Section [IV B 7| summarizes our findings. Where 
applicable, we compare the chirally extrapolated results to experiment. 



1. Isovector Dirac form factor F^(Q 2 ) 

This section covers the isovector Dirac form factor, F^(Q 2 ). Phenomenologically, this form factor is commonly fit 
using a dipole form at fixed pion mass. We will thus first attempt to fit Fi(Q 2 ) using the dipole form and study the 
stability of this fit as a function of the Q 2 range. Next, we will perform chiral fits using the small scale expansion 
(SSE), Ref [53], which includes explicit A (1232) degrees of freedom [SSJ[55]. We will first compare the expansion 
applied to the Dirac radii, (r 2 ), obtained from the previous dipole fits. We will then study the covariant baryon chiral 
perturbation theory expansion (BChPT) for the same quantity, see Ref. [§1 167j. 

Finally, we will present SSE fits to the simultaneous Q 2 and dependence of our lattice data. The latter method 
has the strong advantage that no reliance on the applicability of the dipole form is assumed. For these fits we apply 
the super jackknife and error correlation matrix methods discussed in Sec. |III C| Thus, we believe that this fit strategy 
is superior to the ones previously employed. 

a. Dipole fits to isovector F± (Q 2 ) In this section we discuss the Q 2 -dependence of the form factors at fixed values 
of the pion mass, m T . The function we will use throughout this section is the dipole formula, 

FV(Q 2 )=A Q /(l + Q 2 /M 2 ) 2 , (37) 



with Aq fixing the overall normalization and Md being the dipole mass. From Eqs. (37 1 and (36) it is immediately 
obvious that the dipole mass is related to the Dirac radius of Fi(Q 2 ) via 



In order to verify if the functional form indeed allows for a meaningful application of the dipole formula, we have 
performed a series of fits in which we varied the fit interval [Q 2 in , Q 2 a J and listed the variation of the fit parameters 



We have restricted ourselves to the 28 3 lattice with pion mass — 356MeV. Results are summarized in Tab. [VI 
The table shows the fit interval used, the resulting value of x 2 /dof (degrees of freedom), the normalization Aq - 
which must be equal to one within precision due to the conservation of the vector current — and the dipole mass, 
Md, together with the resulting Dirac radius, (r 2 ). All error estimates have been obtained by app lying the Jackknife 



method to the minimization of the \ 2 including the error correlation matrix, as discussed in Sec. Ill C The table is 
divided into three blocks — first, the large-Q 2 cut-off is varied, next the small-Q 2 cut-off is varied and finally the 
fit-interval window is moved along the available data set. Note that when one leaves out the small Q 2 values, the 
data point Fi(0) is no longer included in the fit interval and Aq can vary more. 

The overall conclusion is that Ao is always compatible with one within error bars and all results for Md are consistent 
over the entire table. The former is an important internal consistency check and the latter allows us to conclude that 
the dipole function is indeed an excellent description of the F^(Q 2 ) form factor over the entire range of available Q 2 
values. 
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[QL^QLJ fGeV 2 l 

[^mm) ^6 max J [ J 


x 2 / dof 


Aq 


M d [GeV] 


(rf) [fm 2 l 

\ -L / L J 


[0, 1.51 


1.22 


1.0004(20) 


1.299(21) 


0.2770(92) 


[0,0.5] 


1.09 


9.9931(21) 


1.298(23) 


0.2774(97) 


[0, 0.41 


1.27 


9.9956(22) 


1.293(23) 


0.2797(99) 


[0 31 


1.52 


9.9932(22) 


1.287(22) 


0.2820(98) 


[0,0.2] 


2.43 


1.0006(23) 


1.276(24) 


0.2870(11) 


[0, 1.5] 


1.22 


1.0004(20) 


1.299(21) 


0.2770(92) 


[0.1,1.5] 


1.23 


0.9965(42) 


1.305(22) 


0.2744(93) 


[0.2,1.5] 


0.72 


0.9780(80) 


1.338(25) 


0.2611(97) 


[0.3,1.5] 


0.90 


0.9765(148) 


1.339(31) 


0.2605(120) 


[0,4,1.5] 


1.21 


0.9445(354) 


1.384(59) 


0.2440(207) 


[0.5,1.5] 


1.41 


0.9469(692) 


1.378(105) 


0.2462(376) 


[0.3,0.5] 


0.65 


1.0098(499) 


1.288(85) 


0.2819(370) 


[0.2,0.4] 


0.77 


0.9749(117) 


1.351(38) 


0.2558(146) 


[0.1,0.3] 


1.59 


0.9948(47) 


1.297(24) 


0.2778(104) 



TABLE VI: Dipole fits to isovector F? (Q 2 ) with varying fit intervals on the 28 3 lattice with m T = 356MeV. 



After performing similar fits at all available pion masses, we obtained the numbers compiled in Tab. VII We have 
taken all available Q 2 values for each fit at fixed pion mass. Again, we have performed a combined error analysis with 



error correlation matrix and jackknife. We find that the resulting values for the Dirac radii, 



are systematically 



m v [MeV] 


X 2 /dof 


A 


M d [GeV] 


(r?) [fm 2 ] 


293 


0.89 


0.9983(73) 


1.307(35) 


0.2734(147) 


356 on 28 3 


1.22 


1.0004(20) 


1.299(21) 


0.2770(92) 


356 on 20 3 


1.95 


0.9999(15) 


1.382(20) 


0.2447(70) 


495 


1.57 


0.9994(9) 


1.3829(12) 


0.2444(42) 


597 


3.76 


0.9998(5) 


1.4144(8) 


0.2336(27) 



TABLE VII: Dipole fits to isovector F?(Q 2 ) for all data sets. 



smaller than the experimental value, 



0.637(12)fm 2 from Ref. [55]. We discuss possible resolutions of this 



discrepancy in the following by discussing chiral fits to our lattice data. 



b. SSE fits to isovector 



As we have seen, the Dirac radius — and, consequently, the size of the nucleon 



- is lower than experiment at the pion masses we use. Hence, for a meaningful comparison the Dirac radii (rf) 
obtained in Tab. VII need to be extrapolated as a function of the pion mass, m v . In this chapter we perform the 



chiral extrapolation using the small-scale expansion (SSE). The pion mass dependence to next-to-leading order (NLO) 
is given by: 



1 



(4^)2 
125[ 



1 



(A) 



(105 



26 



-»*(x) 




(39) 



This expansion has a logarithmic divergence at — ¥ 0. Since the proton and the neutron are linear combinations 
of the isovector and isoscalar operators, the radius (rf) of either the proton or the neutron will similarly diverge in 
the chiral limit. In our fits we fix some of the parameters involved in the chiral expressions, see Tab. | VlTT| The axial 
coupling gA in the chiral limit has been set equal to 1.2 according to Refs. [H [49j [69]. The chiral limit values of 
f„ and the nucleon mass have been determined in Ref. [70] and Refs. [2TJ [7TJ [72] , res P ec ti ver y- The delta-nucleon 
mass splitting A is taken equal to its physical value from the position of the delta resonance pole. Without any loss 
of generality, we set the regularization scale A equal to lGeV. As input for our fits of (rf) we take the Dirac radii 
obtained from the dipole fits in Tab. VII We perform a two parameter fit with the counterterm Bl (X = lGeV) and 
the coupling ca- We vary the upper cut-off in m v and collect the resulting fits in Tab. |LX| The quality of the fits 



is not so good; we will see in Fig. U0\ that the curvature of the SSE curve is stronger than that of the lattice data. 
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Parameter Value 

9A 1.2 

ca 1.5 or free parameter 
U [MeV] 86.2 

mjv [GeV] 0.8900 

A [GeV] 0.2711 

9i 2.5 



TABLE VIII: Input parameters used for the chiral expansion. 



Furthermore, the value of ca tends to be larger than 1.5. One interpretation is that the range of validity of the SSE 
does not extend to our lattice data. On the other hand, if the smallest data point suffers from finite-size effects, it is 
still possible that the SSE is consistent with lattice data on a very large volume. 



m-n max [MeV] 


X 2 /dof 


CA 


Bfo(lGeV) 


500 


8.0 


1.951(36) 


1.713(78) 


600 


7.2 


1.873(20) 


1.557(48) 



TABLE IX: NLO SSE fits to the isovector Dirac radii, (rf), obtained from dipole fits. 



c. BChPT fits to isovector 



In the covariant baryon chiral perturbation theory (BChPT) scheme we use, the 



A (1232) degrees of freedom are not explicitly included while recoil corrections to the non-relativistic Heavy Baryon 
results which correspond to kinetic insertions in the nucleon propagator are systematically resummed [73] . The chiral 
expansion of the (rf) Dirac radius, is given by |67) : 



(r()=B cl + ((r())^ + ((r())W+0{p b ) 1 
The notation ((•)) ( - M - ) denotes the contribution of a quantity to the Mth order in the expansion in p M . 



(40) 

The specific 



expressions are 



B c i 

-2\\(3) 



124(A), 
1 

16tt 2 /2m 



15<7>2 M 



2(5g 2 A + l)A4l g 



A 



M 



N 



2 „2 /i r„,2 



B M N ) 



+ 



g A m * 



«r?» 



(4) _ 



3c 6 g A ml 



fl5ml - 74mlM 2 



A? 



70M 



arccos 



2M N ' 



16^/2^/4^4^2" 



mt 



m 7r (m 2 . — 3Mg) arccos 



2Mn 



4M 2 - ml 



M 2 + (Af 2 -m 2 )log 



(41) 



The expression up to order 0(p A ) introduces two fit parameters, d r & {\) and cq. We set the scale of dimensional 
regularization A equal to the value of the nucleon mass in the chiral limit, M . We point out that when only the 
third order in expansion (W0| is considered, one should replace Mn by Mq in Eqs. (41 1. In our fits, however, we only 



consider the full expansion at order 0(p 4 ). Consequently, we are always using a pion-mass dependent form, Ma^to^), 
for Mat. Similar to Ref. [9] we use the expansion from Ref. [74] to model the functional form Mjv(m 7r ): 

,2 



+ 



3g 2 A ml 



3m* 
"128^7! 

+4e;(AK 



-4 



A 



(>2 



ml/M 2 
M 

87T 2 /2M 2 ° g 



M 2 

A 



■4ci 



M 3 



arccos 



™7T 

2M 



8c! + c 2 + 4c 3 log 



M 



(42) 
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with ci, C2, C3, and eJ(A) being parameters that need to be fixed. In order to fix these constants we follow Ref. [5] 
again and pick the parameters C2, C3, and C4 from the literature, cf. Refs. [72l I75H75] . For M we could adopt the 
chiral value tum listed in Tab. |VIII| which we have used for the SSE. However, since the functional dependence on 
M is highly complicated and all expressions can be sensitive to small changes in M , we decided to consider it a 
free parameter and determine both Mq and the remaining parameters ci and eJ(A) from a fit to the data points 
listed in Tab. [TT] with pion masses m v < 500MeV. The physical nucleon mass is included in the fit. The resulting 
X 2 /dof=2.93 is acceptable and the parameters are in agreement with those reported in [9\. If we include all nucleon 
masses up to m ff < 600MeV, the x 2 /dof=4.15 increases slightly. Since we do not know the exact range of validity, 
we adopt the more conservative choice < 500MeV. When the experimental point is not included, the uncertainty 
of Mo increases dramatically and becomes inconsistent with Ref. [9]. The results are summarized in Tab. [X] As an 
alternative to Eq. (42) we may use an expression from Ref. [21]. It turns out, however, that the influence of different 
expansions of MN(m 7r ) on the radii is quite small. 



M [MeV] 


e{(\ = lGeV) [GeV" 3 ] 


ci [GeV" 1 ] 


c 2 [GeV" 1 ] 


c 3 [GeV" 1 ] 


c 4 [GeV" 1 ] 


876.6(17) 


1.27(12) 


-0.983(22) 


3.2 


-3.4 


3.5 



TABLE X: Low energy constants involved in the chiral expansion of the nucleon mass up to NNLO. Mo, ei(lGeV), and ci are 
determined from a fit to lattice data and experiment. 



There is one more subtlety that prevents a straightforward BChPT fit in the way it was possible for the SSE: The 
parameter cq in Eq. (41 1 appears only at £>(p ). The physical meaning of this low-energy constant is that of the 



chiral limit of the isovector anomalous magnetic moment. In order to gain statistics and better constrain our free 
parameters, we follow the prescriptions devised in Ref. [5] and perform a simultaneous fit to (r 2 ), (r?,), and k v in the 
BChPT scheme. This is done later in the BChPT part of Sec. |IVB 2| and the results are summarized in Tab. |XIX| 
Note that we restrict the joint fit to pion masses m T < 400MeV. The reason is that both the (r 2 ) and the (r|) have 
a stronger curvature at larger pion masses than our lattice data and the x 2 /dof would be unacceptably large when 
pion masses beyond 400MeV are included. We cannot exclude that the range of validity of the BChPT is smaller 
than that of the SSE. Figure [lO] displays the results from the fits in Tabs. |lX| and |XIX| with < 500MeV (for the 
SSE) and < 400MeV (for the BChPT) graphically. The experimental value in the graph is taken from Ref. 68J. 
Both curves are quite similar near the physical point. However, the lattice data is still quite flat for the available pion 
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FIG. 10: Chiral fits to isovector Dirac radii, (rf), for < 500MeV (SSE) and m„ < 400MeV (BChPT). 

masses. Our lattice results do not give any indication of the divergence we should see approaching the chiral limit. 
We also remark that we observe a finite- volume dependence at = 356MeV of almost 12% if we perform dipole fits 



to all available Q 2 , as discussed in Sec. HIE This effect could potentially increase rapidly at lower pion masses, and 
therefore the smallest ~ 293MeV data point could be too low. Thus, it is not excluded that the fit works well 
below 400MeV if one had data at larger volume. At this point we observe that the effect goes in the right direction 
and is qualitatively consistent. For a conclusive statement, we still require more data at smaller pion masses and at 
larger volumes. 
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d. Simultaneous expansion fit to isovector Fi(Q 2 , m w ) The dipole fits to Fi(Q 2 ) discussed previously indeed 
provide an excellent description of the data over the entire range of available data. However, applying dipole fits 
suffers from the disadvantage that there is little fundamental justification for their usage. They are simply employed 
because they appear to work well. As an alternative strategy, we consider now the simultaneous chiral expansions in 
both m n and Q 2 . These expansions are only expected to hold for small values of Q 2 and m n , but we can apply them 
to our entire data set, i.e. by combining data from several different ensembles at different m n values in the same chiral 
fit without model-dependent assumptions on the functional form. 

In the following, we present results from the application of the expression from Ref. [53] to the isovector form factor 
F?(Q 2 ,m n ): 



Fl{Q 2 ,m^) = 1 



1 

dx 



Q 2 





+ f dx 
Jo 



^4-^-2^)+ log( T 



1 ®A 2 c 2 A + ml [3g A + l- 



fciQMl-) A1 ° Si?,,i ' 



V A 2 — rh 2 



32 

dx — c 2 A A y/ A 2 — m 2 log i?(m w ) — \/A 2 — m 2 log R(r 
3 



log 



(43) 



where we use the function R(z) defined in Eq. (30) and introduce 

%l - Q 2 x(l - x) . 



rh 2 



(44) 



Note that the expansion in Eq. (43) is finite at Q 2 = 0, i.e. F^(Q 2 — 0, m OT ) = 1 for any value of including zero, 
but its derivative w.r.t. Q 2 at the origin will diverge logarithmically as — > 0, cf. Eq. (39 1 . Again, the fit parameters 
are i?[ (lGeV) and ca- The other parameters are fixed at their chiral values, see Tab. 



VIII 



As we will also point out in Sec. |IV B 3[ our comparatively large lattice volume puts us at a disadvantage when 
studying the ratio F 2 (Q 2 ) / F 1 (Q 2 ) for large values of Q 2 > lGeV 2 . However, when using chiral expansions, this choice 
turns out to work to our advantage since we have sufficiently many data points at lower values of Q 2 to meaningfully 
apply the chiral expansion directly to the Q 2 -dependence. In the chiral fits we will certainly need to apply cuts in 
Q 2 . Hence, every fit will have fewer data points contributing to it than in the dipole case. Thus, in the end the 
uncertainty may very well turn out to be larger. However, this is offset by the advantage that we do not make any 
phenomenological assumptions on the Q 2 -dependence and the results can truly be considered "first principle" results. 

As the first step, we fix the interval in Q 2 to be Q 2 £ [0,0.5] GeV 2 . This choice is a reasonable guess based on 



the discussion in Ref. 



With this interval, we vary the cut representing the upper value of the pion mass, 



The resulting values for \ /dof, the fit parameters and the extrapolated Dirac radii, (r 2 ) at the physical pion mass, 



are listed in Tab. XI From this table we conclude that it is necessary to apply a rather conservative cut and restrict 



m, max [MeV] 


X 2 /dof 


5io(lGeV) 


CA 


(rl){mt y l [fm 2 ] 


300 


0.59 


0.60(39) 


1.15(33) 


0.623(68) 


400 


1.86 


0.35(11) 


1.009(97) 


0.686(21) 


500 


12.74 


1.195(46) 


1.623(23) 


0.5355(98) 


600 


30.74 


1.181(29) 


1.631(13) 


0.5446(65) 



TABLE XI: Fits to isovector Ff(Q 2 ,m w ) at fixed interval of Q 2 = [0,0.5] GeV 2 with different pion mass cuts. The experimental 
form factors are not included in the fit. 



ourselves to m x < 400MeV. The value of x 2 /dof becomes extremely large beyond this point which implies that the 
function fails to properly describe the m, dependence of the data. This fact is not surprising and has been reported 
also in [5SJ [TSJ [5D] . We have also observed it above when applying the SSE and BChPT formulae to the radii obtained 
from dipole fits. In all those cases it became evident that the pion mass dependence of the lattice data is weaker than 
NLO SSE demands. Although the values of the fit parameters in Tab. XI stabilize when ensembles at > 400MeV 
are included, the resulting parameters are incompatible with those obtained in the region < 400MeV. Next, we 
need to ascertain that the cut in Q 2 can be justified. We vary the upper end of the fit interval in Q 2 by keeping 
the cut for < 400MeV in place. The results are shown in Tab. |XII| It turns out that the fit quality only mildly 
depends on the cut we apply in Q 2 . This situation is quite different from the pion mass dependence we encountered 
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Q 2 max [GeV 2 ] 


xVdof 


B{ (lGeV) 


CA 


(r?)(m* s ) [fm 2 ] 


0.7 


1.84 


0.483(71) 


1.127(54) 


0.661(14) 


0.5 


1.86 


0.35(11) 


1.009(97) 


0.686(21) 


0.3 


1.82 


0.60(21) 


1.241(158) 


0.645(37) 



TABLE XII: Fits to isovector Fi(Q 2 ,m^) at varying intervals of Q 2 with fixed pion mass cut, m n < 400MeV. 



previously. The SSE expression Eq. (43) provides an excellent fit to the Q 2 dependence for essentially the entire range 



of data points, but fails to describe the pion mass dependence for all except the smallest masses. 

We now discuss the resulting curves for the chiral expansion with the cuts Q 2 < 0.5GeV 2 and < 400MeV. 
Figu re \TT\ shows the data for the ensemble at = 293MeV and the curve based on the best fit parameters listed in 
Tab. XII applied to the same pion mass. Regarding the Q 2 -dependence in Fig. 11 we observe the surprising feature, 
that the resulting curve appears to fit the data very well over a very large range of Q 2 values. However, we do not 
believe this to be of physical significance since there is no reason to believe in the validity of the chiral expansion 
at Q 2 values as large as lGeV 2 ! Hence, we consider this feature to be merely accidental and the upper cut in Q 2 
necessary on theoretical grounds. 




FIG. 11: SSE best fit to isovector form factor F"(Q 2 , m w = 293MeV) with lattice data for that ensemble. 



Figure 12 shows the resulting chiral extrapolation of the Dirac radii as a function of the pion mass, m„. For 
illustration purposes, we have included the radii obtained from the dipole fits, cf. Tab. |VII[ in the plot. These data 
points have no influence on the curve presented and just serve as a comparison of the two fitting methods. The red 



star in the plot shows the experimental value taken from Ref. [57] • When studying the resulting Dirac radii in Fig. 12 
we find that the fit slightly overshoots the experiment and that the lattice data may very well account for the physical 
value of (r 2 ). Finally, we point out that the fitted parameters ca and -B[ (lGeV) from the simultaneous strategy are 
different from those obtained with the previous strategy, cf. Tab. |IX| 

From the resulting % 2 /dof we do not see that the SSE at NLO fails to describe the functional form of Fi(Q 2 ,m v ) if 
we consider the parameter region Q 2 < 0.5GeV 2 and < 400MeV. The simultaneous fit yields larger uncertainties 
than an SSE or BChPT fit to Dirac radii from dipole fits and prefers a smaller value of ca- In general a combined fit 
has the advantage that no phcnomcnological assumption on the functional behavior w.r.t. Q 2 is needed. Therefore 
this approach should become the method of choice as sufficiently small and Q 2 are reached. 



2. Isovector Fault form factor F£ (Q 2 ) 

Since according to the definition Eq. ([8]), F%(Q 2 ) involves a spin-flip, this form factor may be better described with 
an additional suppression of Q 2 . On the other hand, as we will discuss in the section on asymptotic scaling, Sec. IV B 3 



the ratio of form factors does not follow the quark-counting rules one would expect from the leading perturbative 
expansion. Hence, one could use either the dipole expression already employed for the form factor Fi(Q 2 ), Eq. (37), 
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FIG. 12: Isovector Dirac radii (rf) with best fit using kinematic cuts Q 2 < 0.5GeV 2 and m„- < 400MeV. The data points shown 
are dipole radii and experiment, which serve to compare this method with the previous one. They have not been included in 
the fit. 



or a tripole expression via 



2\3 



F2{Q')=A /{l + Q'/M() 



with the tripole mass, M t . The tripole mass is related to the Pauli radius via 

™ ~ M 2 • 



(45) 



(46) 



In the following we will first fit F%(Q 2 ) using the dipole form at fixed pion mass and study the stability of this fit as 
a function of cut-offs. Next, we also apply the tripole form to F%(Q 2 ) and study whether the data favors one of the 
two forms. 

Again, we first consider the lattice results at the large volume, 28 3 , at pion mass — 356MeV. At this point 
as well as in the remainder of this section we determine the unc ertain ty by applying the jackknifc method to the 
minimization of \ 2 including the error correlation matrix, cf. Sec. Ill C When comparing dipole and tripole fits for 
the entire range of available data, we find the results listed in Tab. XIII| The fits are applied to the entire range of 
available Q 2 values. They are shown graphically in Fig. 13 



Dipole Tripole 



X /dof 1.04 1.15 
A 3.107(71) 3.044(68) 
M d/t [GeV] 1.067(17) 1.374(20) 
(rf) [fm 2 ] 0.411(13) 0.371(11) 



TABLE XIII: Comparison of dipole and tripole fit for isovector F£(Q 2 ) on the 28 3 lattice at = 356MeV. 



Both results agree within error bars for Ag. The resulting x 2 /dof is almost identical, but the resulting Pauli radii 
disagree by two sigma. Based on the quality of the fits it is not possible to favor either choice since it appears that 
the distinction only becomes important outside of the range of the available data. 

In order to study if both fit strategies remain stable over the entire range of data we pe rform the same variation of 
the fitting interval as previously for Fi(Q 2 ). The results for the dipole fit are shown in Tab. 
tripole fit results in Tab. |XV| 



XIV 



and the corresponding 



When varying the upper cut-off it is apparent that the variation is minimal and the resulting parameters only 
weakly depend on the cut-off. The uncertainty increases as expected, but the effect is small. On the other hand, when 
varying the lower cut-off the error bar increases notably. However, the central values remain stable and the data is 
still well described by the dipole fit over the entire range. When shifting the fit interval, the inclusion of data points 
at smaller values of Q 2 improves the quality of the fit and reduces the error bars notably. In the case of the tripole 
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FIG. 13: Comparison of dipole and tripole fit for isovector F%(Q 2 ) on the 28 3 lattice at = 356MeV. 



[QL„,QLx] [GeV] 2 


X 2 /dof 


A 


M d [GeV] 


(rl) [fm 2 ] 


[0,1.5] 


1.04 


3.107(71) 


1.067(17) 


0.411(13) 


[0,0.5] 


1.15 


3.111(74) 


1.067(20) 


0.411(16) 


[0,0.4] 


1.19 


3.066(78) 


1.095(30) 


0.389(21) 


[0,0.3] 


1.54 


3.065(82) 


1.086(34) 


0.396(25) 


[0,1.5] 


1.04 


3.107(71) 


1.067(17) 


0.411(13) 


[0.2,1.5] 


1.16 


3.113(75) 


1.064(18) 


0.413(14) 


[0.3, 1.5] 


1.54 


3.145(126) 


1.052(33) 


0.422(27) 


[0.4,1.5] 


1.57 


2.824(239) 


1.139(76) 


0.360(48) 


[0.5, 1.5] 


2.05 


2.843(312) 


1.130(96) 


0.366(63) 


[0.3,0.5] 


1.83 


3.866(444) 


0.893(75) 


0.586(98) 


[0.2,0.4] 


1.59 


3.054(92) 


1.098(37) 


0.387(26) 


[0.1,0.3] 


1.54 


3.065(82) 


1.086(34) 


0.396(25) 



TABLE XIV: Dipole fits to F£(Q 2 ) with varying fit intervals. 



fit, the value of A is systematically lower, but still within error bars. We conclude that F^iQ 2 ) is well described 
by either functional form over the entire kinematic range, but the strongest influence on fixing the parameters of the 
fit comes from the region of smaller Q 2 values, in particular from Q 2 < 0.3GcV 2 . Compared to the fits to F^(Q 2 ) 
we thus find a qualitatively similar picture, although the sensitivity of the parameters to the inclusion of the data at 
small values of Q 2 is larger. 

Finally, we perform a series of tripole fits to our entire data set at all pion masses. Results are summarized in 
Tab. XVI We again observe finite-size effects for the Pauli radius (r 2 ). On top of this, the Pauli radius increases as 
the pion mass decreases on the lattices with fixed physical volume. Hence, we have reason to expect that finite-size 
effects may be non-negligible for the smallest value of m^. 

a. SSE chiral fits to isovector (r|) and k v Similar to what we did for the isovector Dirac form factor, we now fit 
the mean squared Pauli radius (r|) and the isovector anomalous magnetic moment, K v , using chiral formulae from 
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[Qmta.QmaJ [GeV] 2 


x 2 / dof 


A 


M t [GeV] 


{rt) [fm 2 l 

\ * 1 L J 


[0, 1.5] 


1.15 


3.044(68) 


1.374(20) 


0.371(11) 


fO 51 


1.15 


3.072(72) 


1.358(24) 


0.380(14) 


[0,0.4] 


1.25 


3.034(76) 


1.386(36) 


0.365(19) 


[0, 0.3] 


1.53 


3.046(80) 


1.364(40) 


0.377(22) 


[0,1.5] 


1.15 


3.044(68) 


1.374(20) 


0.371(11) 


[0.2,1.5] 


1.24 


3.029(71) 


1.377(21) 


0.370(11) 


[0.3,1.5] 


1.60 


2.998(108) 


1.384(38) 


0.366(20) 


[0.4,1.5] 


1.57 


2.709(204) 


1.490(88) 


0.316(37) 


[0.5,1.5] 


2.04 


2.714(261) 


1.483(11) 


0.319(47) 


[0.3,0.5] 


1.84 


3.643(361) 


1.183(85) 


0.501(72) 


[0.2,0.4] 


1.59 


3.004(86) 


1.400(43) 


0.357(22) 


[0.1,0.3] 


1.53 


3.046(80) 


1.364(40) 


0.377(22) 



TABLE XV: Tripole fits to F 2 " (Q 2 ) with varying fit intervals. 



rriir [MeV] 


X 2 /dof 


A 


M t [GeV] 


(ra> [fm] 


293 


1.31 


2.896(162) 


1.389(53) 


0.363(28) 


356 on 28 3 


1.15 


3.044(68) 


1.374(20) 


0.371(11) 


356 on 20 3 


2.09 


2.958(76) 


1.436(25) 


0.340(12) 


495 


1.39 


3.210(44) 


1.482(14) 


0.319(6) 


597 


1.93 


3.402(27) 


1.529(94) 


0.300(4) 



TABLE XVI: Tripole fits to isovector F%(Q 2 ) for all data sets. 



the SSE expansion at NLO. In this case, the expressions as functions of m w are [66), 181]: 



(rl) 



K) 



g 2 A m N 



c\m N 



f%K v (m n )Trmv $f%K v {m n )Tp^/~& - ml 




9A nl ^ m N . 
K„ : — ^ 1- 



2c^tAmjv 



^ 2 f 2 



^ log R(m v ) + log 



BC1 ^/^^ 2, ^c A cv9Am N ml /2A 
-8E[ (X)m N m„ + log I — 



4cAC V gAm N ml 



8c A c v g A A 2 m N 




2 \ 3/2 



log R(m n ) 



277T/2A 

3ml 



24m N 
Kyi^m^) 



2A 2 



(47) 



where R(m) was defined in Eq. (30 1. The pion mass dependence of the product K„(m w )(r|) contains a single unknown 
parameter, B C 2, whose effect is similar to the counterterm B[ (A) in (r\). In addition to the expression above, we need 



to correct for the fact that Eq. ( 47 ) assumes the nucleon mass to be constant when computing the m^ dependence of 
K v . This assumption is certainly not justified, so we need to correct the expression by working with the normalized 
anomalous magnetic moment, K" orm , as defined in Eq. (35). The SSE pion mass dependence for k v involves three 
additional low-energy constants, cy, and El(X). Since K v (m n ) is also part of the expression for (r|) and is 
always obtained from the same tripole fit to F^Q 2 ) it makes sense to fit k v separately at first and then perform a 
single-parameter fit to n v (r 2 ). In this way, the best fit stability is guaranteed. 



We employ the same parameters for the expression as listed previously in Tab. VIII In this case we fix c A to the 
value of 1.5. Furthermore, s ince k v has three fit parameters, we can only work with a pion mass cut of m n < 600MeV. 

The subsequent fit to ^(rf) with n v being canceled out in the next step is shown 



The result is shown in Tab 
in Tab 



XVIII 



XVII 



Unfortunately, the resulting \ /dof is not very good, evidently due to the stronger curvature of the 



NLO SSE form. The curves will be plotted below in Fig. nil 
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m w max [MeV] 


X 2 /dof 


K v 


c v [GeV -1 ] 


£T(lGeV) [GeV~ 3 ] 


600 


1.18 


4.68(24) 


-2.84(37) 


-6.08(41) 



TABLE XVII: SSE fit to the isovector anomalous magnetic moment, K£ orm , obtained from tripole fits. 



m,r max [MeV] 


X 2 /dof 


B c2 [GeV" 3 ] 


600 


22.2 


-6.1(12) 10" 2 



TABLE XVIII: SSE fit to the isovector Pauli radius, (rf), obtained from tripole fits. 



b. BChPT fits to isovector (r|) and k v The expansions of the isovector anomalous magnetic moment, k v , and 
the isovector mean squared Pauli radius, (r|), up to the order 0(p A ) are given by: 
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with the individual contributions at each order: 
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(48) 
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(49) 



In the Pauli radius, eg only shows up in the 0(p ) part which is of higher order. Thus, we cannot use the mean squared 
for a determination of c§ in isolation. We therefore settle for a simultaneous fit to all three quantities 



radius 



(r 2 



(r 2 ) and k v . To do this, we insert the expansion from Eq. (|42| for the nucleon mass Mjy 



in Eqs. (|40f, (|41), 

(481 and (49). The parameter c 4 is again taken from Tab. [Xj We take correlations into account by using both the 
error-correlation matrix and the super-jackknife prescriptions as outlined in Sec. IIIC For these fits we also set the 
scale to A = Mq. The resulting fits are stable, however the x 2 /dof is relatively large. The reason is that the BChPT 
curves for both (r 2 ) and (r|) show a stronger curvature than our lattice da ta. If we included data at higher pion 
masses, x 2 /dof would increase even more. The results are summarized in Tab. XIX Both the (poor) quality of the fit 



X 2 /dof 


C6 


dg(Mo) [GeV" 2 ] 


e? 4 (M ) [GeV" 3 ] 


eWMo) [GeV" 3 ] 


10.23 


4.31(11) 


0.924(12) 


1.201(34) 


-0.08(11) 



TABLE XIX: Simultaneous fit of the BChPT expression to the anomalous magnetic moment, k v , and the radii (rf) and {r^}, 
as obtained from tripole fits. 



as well as the parameters agree with Ref. [9]. We finally display the fits with m^ < 600MeV (SSE) and m w < 400McV 
(BChPT) for the normalized isovector anomalous magnetic moment, K" orm , on the left panel of Fig. 14 and the results 



from the fits with m^ < 600MeV (SSE) and m^ < 400 MeV(BChPT) for the Pauli radius, (r 2 ), in the right panel of 
said figure. 

We find that the fit quality for /c" orm is good and the fits to (r 2 ) are similar to those for (r 2 ) — the lattice data 
does not yet display the feature of strong divergence in the chiral limit and tends to be flatter as a function of the 
pion mass than the associated curves from chiral perturbation theory. Again, it is possible that the fits work already 
well below m n = 400MeV for calculations carried out in sufficiently large spatial volumes. 
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FIG. 14: Results of SSE and BChPT fits to the normalized anomalous magnetic moment 
mean squared radius (rl) (right panel). 



(left panel) and to the isovector 



c. Simultaneous fits to F£(Q 1 , m^) For the chiral expansion in F%{Q , m^) we use the simultaneous expansion 
based on SSE in Q 2 and m w from Ref. [Sol IBT] : 



32c 2 A m N A f 1 
' 9(4tt/ w )2 



dx 



i l0S V 4A-' 



VA 2 - rh 2 y/A2 - 
log iJ(m) - ^ log %) 



(50) 



with i?(a;) and m as defined previously in Eq. (30 1 and (44 1 
for c„(ffl,) defined previously in Eq. (47) 
£[(A 



Note that this expansion also includes the function 

~ The 



VIII 



parameters are Cy, E\(\ — lGeV), and K 
but we decided nevertheless not to add the experimental data point i n this fit 



The input parameters are the same as used previously, Tab 
in addition to ca- The number of parameters is larger than for F^{Q 2 



Sec. 



IVB1 



we also vary the cuts in Q and to find an 



Similar to our study for the form factor -Ff(Q 2 , 
acceptable fit range for F£(Q 2 ,m n ). We start out by keeping the cut Q 2 < 0.4GeV 2 in place and varying the cut in 
the pion mass. Table |XX| summarizes our findings. Note, that we cannot use m n < 400MeV since we would have 



insufficient data points to constrain all three parameters in k v , cf. Eq. (47). As seen from the resulting values of x /dof 



m-n max [MeV] 


x 2 /d°f 


c v [GeV" 1 ] 


SIXiGeV) [GeV" 3 ] 


K v 


ca 


500 


1.82 


-4.1(20) 


-4.7(11) 


4.47(50) 


0.852(98) 


600 


1.61 


-3.54(46) 


-4.40(31) 


4.62(14) 


0.851(97) 



TABLE XX: Fits to isovector F£(Q 2 



at fixed interval of Q = [0, 0.4] GeV with different pion mass cuts. 



the quality of the fits is good for the entire range of pion masses. This feature is distinct from the correspondin g case 
of Fi(Q 2 ,m w ). Next, we vary the upper cut in Q 2 and keep < 500MeV. The results are collected in Tab. XXI 
When increasing the number of data points, the x 2 /dof increases and we conclude that the selection Q 2 < 0.4GeV 
is preferred. Apparently, the Q 2 dependence is not described well beyond that point. We also remind the reader that 



we find finite- volume effects to be reduced when only considering Q 2 < 0.4GeV 2 in Sec. HIE Hence, we finally settle 
for the cuts Q 2 — [0, 0.4] GeV 2 and < 500MeV for our final plots. This choice is compatible with our choice for 
P\i.Q j^w) an d it gives both acceptable fits and minimizes the unknown higher-order contributions in the SSE. We 
notice, however, that the parameter ca is systematically lower than in the previous SSE fits and not compatible with 
its estimate from the A — > Nit decay width 0. Graphically, the resulting curve is shown for a single ensemble at 



rn T[ 



293MeV in Fig. 15 The graph makes it apparent that F^iQ 2 ^^) should be cut below Q 2 < 0.4GeV 2 



We finally plot curves of (rf } and 



1G 



The tripole fit data 

points and the experimental result are included in the graph to provide a comparison to the former method. They have 



derived from the above fit as functions of m n in Fig. 
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g 2 max [GeV 2 ] 


x 2 /dof 


c v [GeV" 1 ] 


E{(lGeV) [GcV~ 3 ] 


K v 


CA 


0.3 


1.54 


-6.6(68) 


-3.8(11) 


4.67(51) 


0.40(32) 


0.4 


1.82 


-4.1(20) 


-4.7(11) 


4.47(50) 


0.852(98) 


0.5 


2.03 


-4.0(14) 


-4.57(78) 


4.63(36) 


0.851(42) 


0.6 


2.77 


-4.4(13) 


-5.11(77) 


4.44(36) 


0.941(35) 


0.7 


4.87 


-4.0(11) 


-5.40(74) 


4.40(35) 


1.100(25) 



TABLE XXI: Fits to isovector F£(Q 2 , m n ) with varying cut in Q 2 at fixed pion mass cut, < 500MeV. 
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FIG. 15: Isovector form factor F$ (Q 2 , m w ) lattice data with best fit SSE at = 293MeV. 



not been used in the fit. The qualitative picture is similar to that from the SSE and BChPT fits to the tripole radii and 
magnetic moments. The curves slightly underestimate the results at the physical pion mass. The advantage of this 
scheme is the same as in the simultaneous fits to Fi(Q 2 , m ff ), namely the absence of phenomenological assumptions 
on the functional behavior w.r.t. Q 2 . 




FIG. 16: Isovector normalized anomalous magnetic moment K" orm (left panel) and mean squared radius (r|) (right panel) as 
a function of the pion mass as obtained from an SSE fit to the chiral expansion of F%(Q 2 , m x ). The data points shown are 
tripole radii and experiment, which serve to compare this method to the previous one and to phenomenology. They have not 
been included in the fit. 
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3. Asymptotic scaling 

Although form factors have been studied experimentally for several decades and using perturbative QCD many qual- 
itative and quantitative features have been understood very well, for the proton F^{Q 2 ) a notable discrepancy has been 
found in recent spin-transfer measurements [59 63J. From quark-counting rules one expects the ratio F| (Q 2 ) / 'Ff (Q 2 ) 
to scale proportionally to Q~ 2 which is consistent with experimental measurements using the Rosenbluth method. The 
recent spin-transfer experiments, on the other hand, found a scaling of F-f {Q 2 )/Ff (Q 2 ) proportional to instead. 
The source of the discrepancy is now generally believed to be two- photon exchange processes, see Ref. [54]. On the 
lattice we are in a unique position to study form factors using exactly single-photon exchanges without contamination 
from other processes. This analysis can proceed in a fully model-independent way. 

The downside of the lattice technology is the limitation to rather small virtualities, Q 2 , since the external momenta, 
p' and p in Eq. ([2]), cannot be chosen too large. For larger values of the external momenta the exponential in Eqs. ( 12 ) 



and (131 introduces large fluctuations which quickly deteriorate the signal-to-noise ratio. A quantitative analysis of 
this phenomenon has been given in Ref. [52] ■ Also, controlling the cutoff effects requires \q\ <C n/a. 

Figure 17 shows our results for the ratio QF% (Q 2 )/Fi(Q 2 ) for all available ensembles. The quantities displayed 



show signs of saturation beyond Q 2 — lGeV 2 for the largest pion masses. For the smaller pion masses saturation is 
not achieved conclusively, but we observe evidence of a pion-mass dependence of the ratio. Unfortunately, we do not 
have data at sufficiently high momentum transfer on the 28 3 lattice at m, = 356MeV to reach the scaling region of 
interest. 
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FIG. 17: Data for the ratio of QF%{Q 2 ) / F?{Q 2 ^ 



To conclude, we find that the picture obtained from spin-transfer experiments is qualitatively reproduced in our 
lattice data. The residual pion mass dependence prohibits a more precise quantitative analysis. New techniques are 
probably needed to establish the relative scaling of the form factors with confidence. 



4- Isovector Sachs form factors 



In the previous Sees. 



IVB1 



and 



IVB2 



we have discussed the Dirac and Pauli form factors, F^(Q 2 ) and F%(Q 2 ) 
and found them to obey dipolc forms very well over their entire parameter range. The Sachs form factors, Ge(Q 2 ) 
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and Gm(Q 2 ), are linear combinations thereof and phenomenologically are usually fit using also dipole forms. It is 
therefore important to see whether the lattice data prefers the one or the other dipole fit scheme since the sum of 
dipoles can only be approximately another dipole form. 

To facilitate this study we have subjected both Ge{Q 2 ) a nd Gm (Q 2 ) to th e same tests as previously applied 
to Fi(Q 2 ) and F% (Q 2 ). The results are summarized in Tabs. XXII and XXIII They show fit results obtained by 

respectively. All tables correspond to the case 



and F% (Q 2 ). The results are summarized in Tabs. XX 
applying various cuts in Q 2 to dipole fits to Ge{Q 2 ) and Gm{Q z 
to_ = 356MeV on the 28 3 lattice. 



[QLn,QL»] [GeV 2 ] 


X 2 /dof 


A 


M d [GeV] 


<r|> [fm 2 ] 


[0, 1.5] 


1.00 


0.9992(20) 


1.029(11) 


0.4406(96) 


[0,0.5] 


0.69 


0.9994(22) 


1.035(12) 


0.4362(103) 


[0,0.4] 


0.80 


0.9997(22) 


1.035(12) 


0.4361(104) 


[0,0.3] 


1.23 


0.9999(22) 


1.035(12) 


0.4364(104) 


[0,0.2] 


1.75 


1.0005(23) 


1.041(13) 


0.4311(111) 


[0, 1.5] 


1.00 


0.9992(20) 


1.029(11) 


0.4406(96) 


[0.1,1.5] 


0.93 


1.0048(42) 


1.022(12) 


0.4474(108) 


[0.2, 1.5] 


0.82 


1.0175(108) 


1.008(16) 


0.4600(150) 


[0.3, 1.5] 


0.81 


1.0500(243) 


0.983(23) 


0.4835(222) 


[0.4, 1.5] 


1.23 


1.0636(625) 


0.975(41) 


0.4915(412) 


[0.3,0.5] 


0.39 


1.0295(649) 


1.003(58) 


0.4649(534) 


[0.2,0.4] 


0.40 


1.0061(151) 


1.023(23) 


0.4461(205) 


[0.1,0.3] 


1.06 


1.0064(51) 


1.024(14) 


0.4455(125) 



TABLE XXII: Dipole fits to G E (Q 2 ) with varying fit intervals. 



[QL„,QL X ] [GeV 2 ] 


X 2 /dof 


A 


M d [GeV] 


(t 2 m) [fm 2 ] 


[0,1.5] 


1.07 


4.068(69) 


1.127(15) 


0.368(10) 


[0,0.5] 


1.21 


4.082(73) 


1.126(18) 


0.369(12) 


[0,0.4] 


1.33 


4.041(77) 


1.148(26) 


0.355(16) 


[0,0.3] 


1.79 


4.051(81) 


1.134(29) 


0.363(19) 


[0.2,1.5] 


1.15 


4.064(73) 


1.127(16) 


0.368(11) 


[0.3,1.5] 


1.55 


4.085(121) 


1.120(29) 


0.372(20) 


[0.4,1.5] 


1.51 


3.716(234) 


1.208(66) 


0.320(35) 


[0.3,0.5] 


1.79 


4.837(415) 


0.966(68) 


0.501(70) 


[0.2,0.4] 


1.70 


4.010(915) 


1.158(33) 


0.348(20) 


[0.1,0.3] 


1.79 


4.051(807) 


1.134(29) 


0.363(19) 



TABLE XXIII: Dipole fits to G M (Q 2 ) with varying fit intervals. 



We also compared the dipole vs. the tripole form for Gm(Q 2 ). The results from the fits are shown in Tab. 



Figure 18 compares the results graphically, 
fitting function. 



XXIV 



Similar to the case of the Pauli form factor we do not find a favorite 



Dipole Tripole 
X 2 /dof 1.07 1.28 

A) 4.0682(694) 3.9976(671) 
M d/t [GeV] 1.127(15) 1.444(18) 
(r 2 M ) [fm 2 ] 0.368(10) 0.336(86) 



TABLE XXIV: Comparison between dipole and tripole fit for isovector Gm(Q 2 )- 



Finally we apply t he dipole fi ts to all available data sets and extract the radii, (rg) and (tm)- The results are 
shown in Tabs. XXV and XXVI It turns out that a dipole fit to Ge(Q 2 ) has a slightly smaller x 2 /dof compared to a 



32 



— Dipole fit 




1- 



0.1 0.2 0.3 0.4 0.5 0.6 
(f[Ge\f] 



FIG. 18: Comparison between dipole and tripole fit for isovector Gm{Q 2 )- 



[MeV] 


X 2 /dof 


A 


M d [GeV] 


{rl) [fm 2 ] 


293 


0.77 


1.0008(74) 


1.014(20) 


0.4544(180) 


356 on 28 3 


1.00 


0.9992(20) 


1.029(11) 


0.4406(96) 


356 on 20 3 


1.66 


0.9998(15) 


1.060(11) 


0.4158(86) 


495 


1.43 


0.9999(9) 


1.106(7) 


0.3818(46) 


597 


4.02 


0.9999(5) 


1.145(4) 


0.3562(28) 



TABLE XXV: Dipole fits to G E (Q 2 ) for all data sets. 



fit to Fi(Q 2 ). On the other hand, fitting Gm(Q 2 ) with a dipole or a tripole is no better than fitting F 2 U (Q 2 ). In both 
cases the fits do not exhibit a clear preference for either a dipole or a tripole function. Similar to the fits to F^(Q 2 ) 
and the case of F%(Q 2 ) we observe finite-size effects when comparing the two volumes at 20 3 and 28 3 . However, in 
the case of Ge(Q 2 ) they appear to be less pronounced, even though the x 2 /dof is slightly smaller. 

5. Flavor dependence 

We investigated the flavor dependence of the connected part of the Fi(Q 2 ) form factor by studying the ratio 
of Ff(Q 2 )/Fi(Q 2 ). Experimentally, this ratio is being scrutinized currently, see Ref. Iffl for the website of the 
experiment. Since the forward value of the ratio is trivially determined by the number of quarks, the interesting 
quantity is the slope, s(m n ), w.r.t Q 2 of the ratio. Before comparing to experiment, this quantity needs to be chirally 
extrapolated and we have adopted the form [S3] 







kx - k 2 log ( 


ml \ 




s 


(47^)2 ; 




[MeV] 


X 2 /dof 


A 


M d [GeV] 


(r 2 M ) [fm 2 ] 


293 
356 on 28 3 
356 on 20 3 

495 

597 


0.93 
1.07 
2.11 
1.13 
1.14 


3.93831(1842) 
4.0682(694) 
3.9854(826) 
4.2913(471) 
4.4684(291) 


1.133(44) 
1.127(15) 
1.180(20) 
1.196(11) 
1.236(8) 


0.364(28) 
0.368(10) 
0.336(11) 
0.3266(62) 
0.3060(38) 



TABLE XXVI: Dipole fits to G M (Q 2 ) for all data sets. 
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with two generic parameters, k\ and k^- Figure [T9| shows the lattice calculations of the ratios at the three lightest 
quark masses. The resulting slopes obtained from linear fits to these data are summarized in Tab. XXVII 




cffGeV 2 ] 



FIG. 19: Linear fits to F?/F?(Q 2 ) ratio. 



m* [MeV] 


Slope [GcV~ 2 ] 


495 


-0.132(5) 


356 on 28 3 


-0.136(11) 


356 on 20 3 


-0.159(9) 


293 


-0.148(18) 


m* s = 139 


-0.190(26) 



TABLE XXVII: Slopes at Q 2 = OGeV 2 of F?/F?(Q 2 ) from linear fits. The last line shows the chirally extrapolated value. 



We observe that the lattice data are quite linear for Q 2 < 0.5GeV 2 and approximately linear up to lGeV 2 . Since 
the leading chiral expansion Eq. (51 ) only tells us how to extrapolate the slope to the physical pion mass, we only fit 

to 
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the slope at low Q 2 values in these fits. Figure 20 shows the chiral extrapolation of the slope displayed in Fig 
the physical pion mass, denoted by the vertical dashed line. The leading chiral singularity is logm^, cf. Sec. IV B 1 
which explains the divergence of the curve as the pion mass goes to zero. The error band propagates the statistical 
errors to the physical mass, and yields the final result for the slope. This final number is shown in the last row of 
Tab. |XXVII| The fit parameters are summarized in Tab. |XXVIII| 

The linear behavior of Ff /F" 
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The prediction that can be meaningfully compared to experiment is shown in Fig. 
as a function of Q 2 departing from Q 2 = 0, extrapolated to the physical pion mass is plotted, including the one-sigma 
error band. 
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FIG. 20: Chiral fit to the slope of the ratio Ft/F?(Q 2 = OGeV 2 ). The physical pion mass is marked by the dashed red line. 



x7dof 


1.5 


fci 


0.093(23) GeW 2 


k 2 


0.023(11) GeV" 2 




-0.190(26) GeV" 2 



TABLE XXVIII: Fit parameters to the slope of the ratio F? / (Q 2 ) as a function of the pion mass. 



6. Isoscalar form factors 



We now turn our attention to the isoscalar form factors. We remind the reader that we have only computed the 
connected Wick contractions, and that the disconnected diagrams remain to be calculated in the future. In the 
following we discuss the isoscalar form factors, Fl(Q 2 ) and F^Q 2 ) separately. 

a. Isoscalar form factor Ff(Q 2 ) We first apply the well-known dipole fits again to the isoscalar form factor 
F((Q 2 ). First, we study the stability of the dipole fits to variations i n the fi t interval like we did before. We focus 
again on the ensemble with m n = 356MeV on the 28 3 lattice. Table IXXIXl summarizes our findings. It is evident 
that the dipole provides an excellent, stable fit over the entire available data range. 

Applying dipole fits to all available pion masses yields the results shown in Tab. |XXX| We find that even at the 
lightest available pion mass the obtained radii underestimate the empirical value (rf) — (0.782fm) 2 quoted from 



u. 



0.3- 



0.2k 



. m = 140 MeV 




0.5 1 
d [Ge\f] 



FIG. 21: Linear slope of F?/F?(Q Z ) at the physical pion mass. 
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K?mm,QL*l fGeV 2 l 

L^v mm 7 max j L J 


x 2 / dof 


A 


M d [GeV] 


(rf) [fm 2 ] 

\ -L / L J 


[0, 1.5] 


1.20 


3.0511(27) 


1096.9(62) 


0.3883(27) 


[0,0.5] 


0.77 


3.0518(27) 


1098.6(65) 


0.3871(46) 


[0 41 


0.90 


3.0521(28) 


1099.0(66) 


0.3868(46) 


[0,0.3] 


0.91 


3.0516(29) 


1070.0(67) 


0.3883(47) 


[0,0.2] 


0.76 


3.0520(29) 


1098.1(68) 


0.3875(48) 


[0,1.5] 


1.20 


3.0511(27) 


1096.9(62) 


0.3883(27) 


[0.1,1.5] 


1.18 


3.0587(66) 


1090.3(81) 


0.3931(59) 


[0.2,1.5] 


1.09 


3.0356(14) 


1103.4(109) 


0.3838(76) 


[0.3,1.5] 


0.57 


3.0141(306) 


1107.9(151) 


0.3807(104) 


[0.4, 1.5] 


0.64 


3.0522(799) 


1097.8(248) 


0.3877(175) 


[0.3,0.5] 


0.30 


3.0309(963) 


1102.9(375) 


0.3842(260) 


[0.2,0.4] 


1.04 


3.0430(204) 


1102.2(135) 


0.3853(98) 


[0.1,0.3] 


0.85 


3.0589(73) 


1090.7(88) 


0.3928(63) 



TABLE XXIX: Dipole fits to isoscalar F((Q 2 ) with varying fit intervals. 



m n [MeV] 


x7dof 


A 


M d [GeV] 


<r?> [fm 2 ] 


293 


1.08 


3.0421(90) 


1086.6(129) 


0.3958(94) 


356 on 28 3 


1.20 


3.0511(27) 


1096.9(62) 


0.3883(27) 


356 on 20 3 


2.46 


3.0538(22) 


1107.3(68) 


0.3811(47) 


495 


1.78 


3.0541(15) 


1152.7(41) 


0.3516(25) 


597 


1.27 


3.0550(9) 


1201.3(29) 


0.3238(16) 



TABLE XXX: Dipole fits to isoscalar Ff(Q 2 ) for all data sets. 



Ref. [66 substantially by a factor of two. A similar discrepancy also appeared in the isovector case, Sec. IV B 1 
but it seems to be even more pronounced here. Another noteworthy fact is the inconsistency of the forward value, 
F((Q 2 — 0), with the expected value of three. This is most likely a cutoff effect. Although the disagreement is only 
mild (about 1%), it is outside the relative error bars. We will comment on this below when discussing the chiral 
expansion. 

Reference |66j discusses the chiral expansion in m n and Q 2 . The next-to-leading order SSE expression 



F^Q 2 ,™,) =3-125! 



Q 2 

(47T/.) 



(52) 



has only a linear Q 2 dependence and no m w dependence. As we have pointed out before, in comparison with the 
isovector form factor Fi(Q 2 , m^), the derivative w.r.t. of the isoscalar one does not diverge in the limit — > 0. 
It merely approaches a constant value. Since the Q 2 de pendence only has a linear part, we need to again restrict 
ourselves to the regime of small values of Q 2 . Table XXX suggests a non-trivial m^-dependence. However, we observe 



that the results for (r 2 ) are still constant within statistics when we restrict ourselves to the region < 400MeV. 
This cut is consistent with what we have found before in the isovector case. 

With the restriction of < 400MeV we perform a chiral fit based on Eq. (52 1. However, when covering the 
entire fit interval, Q 2 = [0,0.5] GeV 2 , we find an unacceptably large x 2 /dof=128. We traced the problem back to 
the value at the origin, Ff(Q 2 = 0). Due to the systematic shift with a small relative error we end up with such a 
huge discrepancy. Hence, we restrict ourselves to the fitting interv al, Q 2 > O.OlGeV 2 , which excludes this data point. 
Varying the upper cut-off in Q 2 yields the results tabulated in Tab. XXXI Surprisingly, the radii are even smaller than 
those obtained from the dipole expression. Furthermore, the overall quality of the fits is poor until Q 2 < 0.2GcV 2 is 
imposed, at which point the resulting radius is still about 20% smaller than the one obtained from the dipole form. 
We conclude that applying the NLO expression ( 52 ) does not give new insight into the problem. 

b. Isoscalar form factor F^Q 2 ) The chiral expansion from Ref. 
F£(Q 2 ,m w ), an expansion of the form 



yields for the isoscalar spin-flip form factor, 



F|(Q 2 , m w ) = k" - 8m N mlE 2 



(53) 
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Q 2 max [GeV 2 ] 


X 2 /dof 


Si 


(r 2 ) [fm 2 ] 


0.4 


42.7 


0.3304(23) 


0.2631(18) 


0.3 


17.9 


0.3493(25) 


0.2782(20) 


0.2 


2.1 


0.3878(47) 


0.3089(38) 



TABLE XXXI: Fits to isoscalar F({Q 2 , m n ) at varying intervals of Q 2 with fixed pion mass cut, m n < 400MeV. 



i.e., the leading Q 2 dependence vanishes and the dependence is quadratic. Table XXXII shows the result of a 
fit of Eq. (531 to our lattice data with the cuts Q 2 = [0,0.5] GeV 2 and w, < 400MeV. The experimental value for 
F|(Q 2 = 0, m„- = mf ys ) taken from Ref. [SH] is listed on the last line. It is evident that the data suggests n° s < 0, 



Fit result 


X 2 /dof 


1.34 




-0.36(25) 


E 2 [GeV- 3 ] 


-0.32(29) 


Extrap. Fi(Q 2 = 0,mf ,s ) 


-0.32(21) 


Exp. F|(Q 2 = 0,mP^ s ) 


-0.360586110(14) 



TABLE XXXII: Fit results to isoscalar form factor Fi(Q 2 ) for the interval Q 2 = [0, 0.5] GeV 2 and < 400MeV. The next-to- 
last line is our extrapolated result at the physical pion mass and the last line lists the experimental value for F^iQ 2 = 0, m£ hys ). 



but we can only fix the order of magnitude since the relative error is about 50%. When plotting the data for the 
ensemble at — 293MeV, we find that the data is indeed flat, albeit with a large relative error. Figure [22] shows 
this graph with the fit function from Tab. XXXII for F 2 S (Q 2 , m, = 293MeV), which is constant w.r.t. Q 2 . 



O.U 



in 



-0.1- 



0.2 



-0.3 



it ii 



II" 



0.2 0.4 0.6 0.8 
d [GeV*] 



FIG. 22: Isoscalar form factor F|(Q 2 ) as a function of Q 2 with SSE fit at = 293MeV. 

Note that the parameters are poorly determined, although the fit quality as indicated by % 2 /dof is acceptable. We 
can extract the order of magnitude of both k° and E%, but with large statistical uncertainty. On the other hand, 
we have found the extrapolated value of F£(Q 2 = 0, mf ys ) to be fully compatible with experiment, having both 
the correct order of magnitude and the correct sign. We point out again that the forward value will not receive 
disconnected contributions, thus the value quoted corresponds directly to the experiment. 



7. Summary of electromagnetic form factors 



For the isovector form factors of the vector current we find that the lattice mean squared radii at present pion masses 
are significantly below the empirical ones. This finding is qualitatively compatible with chiral perturbation theory, 
though, since the latter predicts a sharp increase of the mean squared radii as the pion mass approaches the chiral 
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limit. Although our data is consistent with one-loop expressions from both SSE and BChPT, we still require data 
at smaller pion masses and larger volumes to successfully compare with phenomenology and make predictions in the 
chiral limit with negligible systematic errors. For the isoscalar form factors of the vector currents we find our results 
to be in agreement with expectations from chiral expansions — the connected part of F((Q 2 ) has a non- vanishing 
Q 2 -dependence, while the Q 2 -dependence of F% (Q 2 ) was not measurable. The overall magnitude of F£(Q 2 — 0) is 
well in agreement with empirical information, albeit with large uncertainty. 



C. Axial form factors 



The nucleon axial current matrix element is also expressed in terms of two form factors: The axial form factor, 
Ga(Q 2 ), and the induced pseudoscalar form factor, Gp(Q 2 ). They correspond to the generalized form factors A 10 (Q 2 ) 
and B W {Q 2 ) in Eq. j^. For a review of the current experimental and theoretical understanding of the axial structure 
of the nucleon, we refer to Ref. [55]. For more details on the chiral effective field theory expressions that we will use 
in our fits, see Ref. [66]. For lattice results from other groups, see Ref. [80] . 

Like the vector current, the axial current needs to be renormalized. The renormalization coefficient can be computed 
from the conserved axial current, see Ref. [53]. Defining the two-point functions C(T) and L(T) of the conserved and 
local currents via 



C(T+ 1/2) = J2(Mx,T+ 1/2Lt(0,0)>, 

X 

L{T) = 5^4o(f,T)7r(ff,0)), 



(54) 



the axial renormalization constant, Za, can then be computed from 

1 / C{T+ 1/2) + C(T- 1/2) 



Za(T) = 



2L(T) 



2C(T + 1/2) 
L(T)+L(T+1) 



(55) 



For large T> 1, the ratio in Eq. ( |55| ) becomes the axial renormalization constant, Za = liniT->oo Za(T). Based on 
the ensembles listed in Tab. [I] wc have used the ranges T = [12, 29] and T = [35, 52] to obtain our final values of 
the renormalization constants. Table XXXIII summarizes our findings. We thus obtain a quark-mass dependent Za 





Volume 


Za 


0.007/0.050 


20 3 


x 64 


1.0816 


0.010/0.050 


28 3 


x 64 


1.0850 


0.010/0.050 


20 3 


x 64 


1.0849 


0.020/0.050 


20 3 


x 64 


1.0986 


0.030/0.050 


20 3 


x 64 


1.1090 



TABLE XXXIII: Renormalization constant of the axial currents. 



factor. Since all the ensembles considered here are at the same value of the lattice spacing a, the differences between 
the Za values in Tab. |XXXIII| are entirely due to quark mass effects. To leading order these effects are expected to 
be linear in am q [87) . and it is important to keep at least this leading quark-mass dependence to avoid introducing 
0{a) discretization errors in the matrix elements of A$. 



G 



We will now r eview t he isovector axial form factor, Ga{Q 2 ), in Sec. IV C 1 the induced ps eudosca lar form factor, 

p(Q 2 ) 



Sec. 



IVC2 



and the two isoscalar axial form factors A W (Q ) and B W (Q ) in Sec. 



summarizes our findings. 



IVC3 



Section 



IVC4 



Isovector axial form factor Ga (Q 2 



The isovector axial form factor is usually fit by a dipole form. The forward value, qa — Ga(Q 2 — 0) has been 
discus sed in de tail in Sec. IV A The dipole mass and its relevance to experiment has been studied in detail in Ref. [85] . 
Table XXXIV shows results of a series of fits to Ga{Q 2 ) on the 28 3 lattice at = 356MeV. It is evident that the 
dipole form provides a good fit to the lattice data at all available values of the virtuality, Q 2 . It turns out that 
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f<2m,„,QLxl [GeV 2 l 

[%mm) max j |_ J 


X 2 /dof 


9a 


M d [GeV] 


(r 2 A ) [fm 2 l 


[0, 1.5] 


1.70 


1.1245(147) 


1.587(29) 


0.1856(67) 


fO 51 


1.77 


1.1334(158) 


1.545(31) 


0.1956(79) 


[0,0.4] 


2.22 


1.1343(160) 


1.536(34) 


0.1980(89) 


[0,0.3] 


1.71 


1.1331(161) 


1.550(39) 


0.1945(97) 


[0,0.2] 


2.25 


1.1308(163) 


1.508(47) 


0.2056(13) 


[0, 1.5] 


1.70 


1.1245(147) 


1.587(29) 


0.1856(67) 


[0.1,1.5] 


1.79 


1.1247(150) 


1.587(29) 


0.1856(68) 


[0.2, 1.5] 


1.43 


1.1303(168) 


1.584(33) 


0.1862(78) 


[0.3, 1.5] 


0.53 


1.1123(228) 


1.598(50) 


0.1830(115) 


[0.3,0.5] 


0.54 


1.1005(642) 


1.617(17) 


0.1787(369) 


[0.2,0.4] 


1.71 


1.1533(209) 


1.496(49) 


0.2088(136) 


[0.1,0.3] 


2.01 


1.1305(171) 


1.563(49) 


0.1913(120) 



TABLE XXXIV: Dipole fits to isovector G A (Q 2 ) with varying fit intervals. 



the axial mass at m, = 356MeV is about 50% larger than the phenomenological value from neutrino scattering, 
M d = (1.026 ±0.021) GeV, and the one from electro- pro duct ion, M d = (1.069 ± 0.016) GeV, see Ref. [85]. Hence, the 
situation is qualitatively similar to the case of the form factors of the vector current. 

For the other ensembles, we obtain the results from the dipole fits at fixed pion mass summarized in Tab. |XXXV| 
The data has been fit to all available Q 2 . The large offset between lattice and phenomenology exhibited by the dipole 
masses, and thus also by the radii, remains to be explained. 



m n [MeV] 


X 2 /dof 


A 


M d [GeV] 


{r\) [fm 2 ] 


293 


0.80 


1.154(26) 


1.577(56) 


0.1879(134) 


356 on 28 3 


1.70 


1.125(15) 


1.587(29) 


0.1856(67) 


356 on 20 3 


1.30 


1.144(15) 


1.661(33) 


0.1694(68) 


495 


1.25 


1.142(77) 


1.654(16) 


0.1708(34) 


597 


1.05 


1.146(46) 


1.686(11) 


0.1644(21) 



TABLE XXXV: Dipole fits to isovector G A (Q 2 ) for all data sets. 



We also study Ga(Q 2 



using the chiral expansion from Ref. [66 . At this order. 



contains a counterterm which is directly related to the mean squared axial radius, 
is encoded in qa (m^): 



the Q 2 -dependence is linear and 
r 2 A ) . The pion mass dependence 



(56) 



For the expression of gA^n) see Eq. (29 1. In the following we have performed fits of the Q -dependence of Ga at 



fixed pion masses and extracted the mean squared axial radii by taking 



-6 dG A (Q 2 ,m„) 



dQ 2 



(57) 



Q 2 =0 



From the resulting series of fits we attempt to extract a sensible interval for the fi tting range of Q 2 . First we determine 
the upper cut-off in Q 2 on the 28 3 lattice. The results are listed in Tab. XXXVI Since the fit expression is only linear, 
the exclusion of points at larger Q 2 does indeed improve the quality of the fit, although the overall result is not as 
satisfactory as for the dipole fit. Both for gA and for (r 2 A ) we find that the chiral fit lies below the results from the 
dipole fits. Furthermore, the fits in this section possess larger uncertainties and the resulting % 2 /dof is worse. 

Now, we apply the SSE fit simultaneously to the entire available data set. From the previous section we have 
learned that Q 2 = [0, 0.4] GeV 2 is a sensible fitting interval which we keep in the following. When varying the upper 
cut-off in m^, we obtain the fit results in Tab. |XXXVli] 

Graphically, Fig. [23] shows the results of the chiral fit together with a dipole fit and the data set for the 28 3 lattice 
at = 356MeV. The fitting range is Q 2 = [0, 0.4] GeV 2 and m ff < 400MeV for the SSE expression and all Q 2 
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Q 2 max [GeV 2 ] 


xVdof 


Bi 


(ri)(m„ = 356MeV) [fm 2 ] 


1.5 


4.82 


0.725(23) 


0.1269(38) 


0.5 


2.66 


0.848(33) 


0.1485(46) 


0.4 


2.67 


0.901(39) 


0.1577(57) 


0.3 


2.52 


0.949(47) 


0.1660(71) 


0.2 


1.50 


1.094(70) 


0.1914(11) 



TABLE XXXVI: Chiral fit to G A (Q 2 ,m^ = 356MeV) with varying upper cut-off in Q 2 . 



vn-r: max [MeV] 


X 2 /dof 


B z 


{r\)(m v = 293MeV) [fm 2 ] 


300 


0.71 


0.993(130) 


0.1707(193) 


400 


1.73 


0.908(32) 


0.1560(60) 


500 


1.54 


0.891(20) 


0.1531(49) 


600 


1.76 


0.847(12) 


0.1457(42) 



TABLE XXXVII: Chiral fits to G A (Q 2 ) at fixed interval Q 2 = [0, 0.4] GeV 2 with varying upper cut-off in re- 



values for the dipole fit. The plot illustrates that the dipole fit works well and that the SSE NLO result is linear 
in Q 2 and provides a good description of lattice data only for relatively small values of Q 2 . There is a rather weak 




FIG. 23: Comparison of dipole and chiral fit to Ga(Q 2 ) for the 28 3 lattice at m n = 356MeV. 



sensitivity to the pion mass, but we need to point out that we do not at all reproduce the value B 3 =3. 08(27) from 
Ref. [53] determined from the experimental value of the axial radius (r\) — 0.409(12)fm 2 (from electroproduction) at 
the physical pion mass. 

To summarize, we find that Ga{Q 2 ) is described excellently using a dipole-type fit, but the axial radius is substan- 
tially smaller than the experimental one. The SSE expansion simply gives a linear dependence in Q 2 — it can thus 
only be applied to small values of Q 2 , where the functional behavior is approximately linear. At fixed m w , the fit 
quality is acceptable. When performing a combined fit in Q 2 and m w , we still find that the axial radius is substantially 
underestimated; we also detect a small residual pion mass dependence that is not present in the SSE expression at 
the given order since at this order one should take gA{ m ir) = 9a = 1<2, cf. Tab. VIII[ in the denominator of Eq. (57l. 
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2. Isovector pseudoscalar form factor 



This section discusses fit results to the induced pseudo-scalar form factor, Gp(Q 2 , m w ). We will first focus on the 
isovector case. Unlike the other form factors, the chiral expansion of Gp(Q 2 ,m 7r ) includes a pion-pole term [66| [80] : 



G P (Q 2 



4m 2 N 
ml + Q 2 



9a 



2ml B 2 
(4tt/ w ) 2 



(58) 



>. The former parameter, B2, is a correction to the residue 



r A ) is the axial radius encountered above in the axial form factor and induces an 



This NLO expression has two fit parameters, B2, and 
of the pion pole. The parameter (r A 
overall shift. 

As discussed previously, we first perform a series of fits on the 28 3 la ttice at fixed m n — 356MeV and vary the 
upper cut-off in Q 2 . For the input parameters to expression (58 1 see Tab 
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This will give us an understanding of 
how far the chiral expansion can be expected to hold, although we always have to ke ep in mind that any agreement 
for Q 2 larger than 0.5GeV 2 should be considered merely accidental. Table XXXVIII summarizes the results. 



Q 2 max [GeV 2 ] 


X 2 /dof 


B 2 


(r\) [fm 2 ] 


1.5 


1.26 


-0.96(13) 


0.107(8) 


0.5 


1.08 


-1.01(17) 


0.108(14) 


0.4 


1.22 


-0.84(22) 


0.088(21) 


0.3 


0.84 


-0.55(31) 


0.051(35) 



TABLE XXXVIII: Pion pole fits to the isovector pseudoscalar form factor Gp(Q 2 ,m n = 356MeV) on the 28 3 lattice. 



The overall quality of the fit is good and % 2 /dof acceptable. However, we notice that the axial radius, (r A ), is 
smaller than the one obtained using the axial form factor, G A (Q 2 , m^). It appears that the NLO SSE formula is not 
able to connect our data to experiment for that observable. 

With a fixed fit interval Q 2 = [0, 0.5] GeV 2 we also performed a combined fit in Q 2 and m n and varied the upper 
cut-off for the latter. Note that in this case, the location of the pion pole varies and is set to the appropriate value for 
each ensemble under consideration. Table XXXIX summarizes our findings. We find the fit formula applicable with 



max [MeV] 


X 2 /dof 


B 2 


(rl) [fm 2 ] 


300 


1.13 


-0.10(80) 


0.058(46) 


400 


1.72 


-0.94(15) 


0.089(12) 


500 


2.85 


-1.544(87) 


0.113(9) 


600 


3.14 


-1.875(56) 


0.139(7) 



TABLE XXXIX: Chiral fits to G P (Q 2 ) at fixed interval Q 2 = [0,0.5] GeV 2 with varying upper cut-off in m*. 



the cuts Q 2 = [0, 0.5] GeV 2 and < 400MeV. 

Next, we study the location of the pole in Eq. (58). Whereas in the previous fits the measured pion mass was used 
to fix the dependence, we can also treat the pole position as a free pa rame ter. For this study, we choose the fit interval 
Q 2 = [0,0.5] GeV 2 on the 28 3 lattice and obtain the results in Tab. 

about 10%, but within this uncertainty, we find that the location of the pion pole is reproduced. 



It is evident that there is an uncertainty of 

We compare the 



Fit result 
X 2 /dof 0.94 
B 2 -1.98(61) 
(ri) [fm 2 ] 0.192(65) 
Fit result m w [MeV] 417(43) 
Actual ensemble m-n [MeV] 356 



TABLE XL: Result from fit of isovector G P (Q 2 ) to pion-pole form with variable pole mass parameter. 
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— it -pole fit, pole position free 
%-pote fit 
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0.6 



FIG. 24: Comparison of pion-pole fits to isovector Gp(Q 2 ) on the 28 3 lattice with fixed pion pole and with the pion pole as a 
free parameter. 



fits with varying pole po sitio n in Tab. 
GeV 2 graphically in Fig. 
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XL 



and with fixed pole position in Tab. XXXVIII for the interval Q — [0, 0.5] 



Our overall conclusion is that the isovector Gp(Q ) form factor can be well described by the pion pole expression. 
We find that by applying the cuts Q 2 — [0, 0.5] GeV 2 and m n < 400MeV we obtain acceptable results. Allowing the 
pion pole to vary freely we obtain a result compatible with the "true" pion mass, albeit slightly higher. 



3. Isoscalar axial form factors 



When considering the forward case of the isoscalar axial form factors, we recover the connected part of the first 



moment of the spin-dependent isoscalar parton distribution which we discuss thoroughly in Sec. IV D Their Q 2 - 
dependence is not known experimentally, but they have been studied in the framework of chiral perturbation theory. 
Reference [55] finds a counter-term with a linear Q 2 -dependence, and Ref. [5j5] does not list any Q 2 -dependence at the 
order considered. In a previous paper [30], however, the same authors find a counter-term with linear Q 2 -dependence, 
in agreement with Ref. [88]. In the following we adopt the notation from Eq. (|9| when referring to the two form 
factors, i.e. we denote them by A W (Q 2 ) and -Bi (<3 2 ). 

a. Isosca lar axial form factor A±a Since the m^-dependence of the forward matrix element is already covered 
we focus on the Q 2 -dependence at each ensemble. First, we study the case = 356MeV on the 28 3 
As a generic fit formula we use the dipole expression. Figure [25] displays the result of the fit. Notice that 
this figure contains an example for six data points that appear superficially high as discussed previously in Ref. [41) . 
cf. also Sec. |III C| The dipole form actually provides a decent description of the data. We find a dipole mass of 
Aid = 1683(89)MeV with x 2 /dof=1.87. Unlike the isovector form factors of the vector current, we cannot refer to 
previous phcnomcnological fits of experimental data. Our results for the dipole fits to the different ensembles are 
listed in Tab. 



in Sec. IV D 
lattice 



XLI 



The radii thus obtained are denoted by 



rrw [MeV] 


X 2 /dof 


A 


M d [GeV] 


(f?) [fm 2 ] 


293 


1.45 


0.560(31) 


1.626(143) 


0.176(31) 


356 on 28 3 


1.87 


0.568(18) 


1.683(89) 


0.165(17) 


356 on 20 3 


1.58 


0.542(17) 


1.903(97) 


0.129(13) 


495 


1.10 


0.586(9) 


1.711(40) 


0.1596(75) 


597 


0.75 


0.605(6) 


1.691(26) 


0.1635(51) 



TABLE XLI: Dipole fits to isoscalar Am(Q 2 ) for all data sets. 



b. Isoscalar axial form factor B\q For the connected part of the isoscalar form factor i?io(Q 2 ) of the axial current 
we also find a non-trivial Q 2 dependence. In this case, however, we observe that the form factor has quite a large 
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magnitude and a rather strong fall-off. Again, we attempt a dipole-type fit and display the result in Fig. |26| The 
dipole appears to fit the connected part well with x 2 /dof=f .88. However, since the functional form is suggestive of 
the pion pole, cf. Eq. (58 1 in Sec. IV C 3 we have also attempted a fit to that form. If we keep the pion pole fixed 
to the actual pion mass of the sample, m n — 356MeV, we obtain a x 2 /dof= 1.86, whereas leaving it a free parameter 
gives x 2 /dof= 1.97 with a measured position of — 394(59)McV. This result is certainly in agreement with the 
actual pion mass of the underlying ensemble. This indicates, that the connected part of the isoscalar i?io(Q 2 ) form 
factor is indeed compatible with a pion pole form, although the data does not favor it over a dipole form. 




0.2 0.4 0.6 0.8 

<f[Ge\f] 



FIG. 26: Dipole fit to the isoscalar axial form factor Bio(Q 2 ) on the 28 3 lattice at = 356MeV. 



The resulting parameters from dipole fits to all data sets are assembled in Tab. XLH It is evident that there is 



[MeV] 


X 2 /dof 


A 


M d [GeV] 


(ft) [fm 2 ] 


293 


0.65 


13.3(37) 


0.671(95) 


1.04(30) 


356 on 28 3 


1.88 


13.14(85) 


0.668(22) 


1.05(7) 


356 on 20 3 


1.22 


11.4(11) 


0.757(36) 


0.82(8) 


495 


1.15 


10.44(54) 


0.802(23) 


0.73(4) 


597 


0.89 


10.27(34) 


0.840(17) 


0.66(3) 



TABLE XLII: Dipole fits to isoscalar Bio(Q 2 ) for all data sets. 
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a notable pion mass dependence, with the dipole mass being affected strongest. This is in line with the functional 
form in Ref. [HS] , although the strength of the effect is quite surprising. At m T = 293MeV the mean squared radius, 
(ff), is larger than lfm, making it the largest radius of all observables. Further studies of the chiral behavior for this 
observable and the computation of disconnected diagrams are interesting directions for future work to shed further 
light on these observations. 



4- Summary of axial form factors 



The isovector axial form factor Ga(Q 2 , Tn^) at present lattice pion masses systematically is much flatter as a function 
of Q 2 than experiment and this result is not explained by one-loop chiral perturbation theory. This is a qualitative 
and significant quantitative mismatch between lattice and phenomenology that still remains to be understood. The 
corresponding pseudoscalar form factor Gp(Q 2 , m w ) is described very well by a pion-pole form and both the functional 
form and the location of the pion pole are in excellent agreement with theory. Still, the SSE can not explain the 
discrepancy in the axial radii we extract from Ga and Gp. The connected parts of the isoscalar form factors of the 
axial current have a strong Q 2 dependence and we have successfully used dipole forms in each case. 



D. Generalized Form Factors 



In this section we present a survey of our results for the generalized form factors, emphasizing their main qualitative 
features and describing the progress made since our previous publication on the subject [20] . 

We start with the form factors of the twist-two quark bilinear operator qD^^/^yq. This operator is particularly 
important since its forward matrix elements determine the quark momentum fractions; in the isosinglet case, it is one 
of the terms appearing in the energy-momentum tensor. Furthermore its forward matrix elements contribute to the 



nucleon mass |91j and spin decompositions, see Sec. IV E The three form factors A20, B20, C20, which determine its 
matrix elements between two arbitrary one-nucleon states via Eq. ([6]), are displayed in Fig. [27] both for the isovector 
and isosinglet combinations. In the latter case, our computation lacks contributions from the disconnected diagrams. 

Similarly, the form factors A 2 q and B 2 q associated with the corresponding twist-two axial operator qD^j u y^/ 5 q, 
see Q, are displayed in Fig. 28 The forward matrix elements of this operator determine the first moments of the 
polarized structure functions, in other words the polarized momentum fraction. As is clear from Eqs. ([6| and Q, only 
the A20 and A20 form factors can be directly obtained at zero momentum transfer; the others require an extrapolation, 
much as the Pauli form factor. 

The statistical uncertainties on A 2 o, B 2 o, C 2 o at = 356, 495 and 597MeV are reduced by about a factor three 
as compared to our previous publication |20| . Furthermore we have data at one lighter pion mass, 293MeV, with an 
accuracy comparable to our previous data in the 600-700MeV range. Given the expected rapid deterioration of the 
signal-to-noise ratio when — > 0, this represents significant progress. 

Similar remarks apply to the axial GFFs A20 and i?2o- The signal-to-noise ratio of the latter GFF is overall quite 
poor and worsens visibly with the pion mass. While previously we were essentially unable to obtain a signal below 
run = 500MeV, we now dispose of some information down to about m, — 350MeV. 

As in our previous work, we observe a qualitative agreement of the relative magnitudes of the GFFs with predictions 
from large N c counting rules, e.g. 

\A u 2 + d \ ~ Nt » l^2o" d | ~ N c , \B^ d \ ~ N? » \B-+ d \ ~ Nl \C%+ d \ ~ N 2 C » \C£ d \ ~ N c . (59) 

As a matter of fact, we remark that C20 is consistent with zero in the isovector channel, whereas it is clearly non- 
vanishing and negative in the isosinglet channel. This is of interest, since C20 entirely determines the longitudinal 
momentum transfer £ dependence of the functions ff™ =2 (£,£) and £/ n=2 (£,i), see Eq. Q. In the case of B 2 o, it is 
the opposite: B 2 q is very small in the isosinglet channel, but it is positive and quite large (compared to A20) in the 
isovector channel. However, B 2 Q d falls off faster than A 2 Q d \ so that the two form factors are practically equal by 
the time \t\ w lGeV 2 is reached, for pion masses < 356MeV. Their t-dependence is thus qualitatively similar to the 
Q 2 -dependence of the electromagnetic form factors, for which chiral perturbation theory predicts that (r^) 2 ~ logm^, 
and (r^) 2 ~ l/m„, and experimentally, it is indeed the case that {r\) 2 < (r^) 2 - 
In the polarized case, the large- N c counting rules predict 

\A%o d \ ~ N 2 » \Al+ d \ ~ N c , \B%- d \ ~ » |B 2 V d | ~ iV c 3 . (60) 

Our lattice data shows that the isovector GFFs are only marginally larger in magnitude than the corresponding 
isosinglet GFFs. We further observe in both the isovector and isosinglet channel that Bw{t) > A2o(t). Noting that 
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FIG. 27: The unpolarized isovector and isosinglet GFFs A20, B20, 
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FIG. 28: The polarized isovector and isosinglet GFFs A20 and i?20- 
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jjn=2 _ an( j jr;n=2 _ see (Igl j g com p a tible with a dominance of the GPD E over at small —t as 
indicated in a recent HERMES study |92| . 

It is also interesting to compare the size of corresponding polarized and unpolarized form factors. The large- N c 
ordering is indeed observed if one compares A%o d with A^ d , or B^ d with B^ d , but the differences in magnitude 
are not pronounced enough to speak of a hierarchy. 

The momentum transfer dependence of the GFFs Aiq, A20, A30 is compared in Fig. 29 for both the isovector and 
isoscalar channels. For this purpose they all have been rescaled to be equal to unity at t = 0. It is clear that in both 
channels, the form factor flattens when the moment n increases. The effect is observed at all pion masses, and it is 
larger when increasing n from 1 to 2 than when going from n = 2 to n = 3. 

This observation can be made more quantitative by fitting the GFFs with a dipole ansatz for \t\ ranging from zero 
to 0.5GeV 2 . In the isovector channel this leads, via the analogue of Eq. (36), to the squared radii displayed in Fig. 30 



The left panel corresponds to unpolarized GFFs, the right panel to the polarized ones. In the infinite momentum 
frame, the transverse radii correspond to the rms transverse distance (i.e. the impact parameter) of the active parton 
to the center of momentum of the nucleon. At each pion mass, the higher the moment of the structure function, the 
smaller the radius associated with it. This effect had been anticipated [35] and was seen by direct lattice calculation 
for the first time in [2]. Indeed most of the contribution to higher moments comes from the large- a; region; a parton 
carrying by itself most of the nucleon's momentum must be located near the center of mass in the transverse plane. 

On the V = 20 3 lattices, corresponding to physical volumes of (2.5fm) 3 , it is only for the electromagnetic radius, 
n = l, that we find a statistically significant increase when the pion mass is reduced. This corresponds to the idea that 
the growth of the nucleon radius with decreasing quark mass is due to the low- a; partons, while the large- x 'valence' 
partons have a transverse distribution that depends only weakly on the quark mass. However, these qualitative lessons 
are questionable due to the statistically significant finite-size effect observed in several of these radii. The larger lattice 
volume V = 28 3 , corresponding to (3.5fm) 3 , leads to larger radii. The largest finite-size effect is seen in the vector 
and axial- vector form factors, see Sec. |IIIE] but a similar effect is also seen in the A20 form factor. It appears that the 
nucleon is somewhat 'squeezed' by the periodic box in which we study it, leading to an underestimation of the radii. 
This finite-size effect is not statistically significant if one simply extracts an effective square radius from the smallest 
available momentum transfer q — 2tt/L, as we already pointed out in [S]. The finite-size effect affects the intermediate 
values of —t most strongly, as discussed in Sec. |III E| We believe that this apparently i-dependent finite-size effect on 
form factors is interesting and deserves further investigation. 

Regarding Fig. [30]we also note that the n = 2 and 3 polarized radii at =356MeV are statistically compatible 
with the corresponding unpolarized ones. This is in contrast with the n = 1 case, where the Dirac radius is significantly 
larger than the axial vector radius. In the chiral limit, the Dirac radius diverges logarithmically, while the axial vector 



radius remains finite, see Sees. IV B and IV C Therefore, the former must become larger than the latter below a 
certain pion mass. What we see on Fig. |30| is that this hierarchy survives up to very large pion masses, far beyond 
the range of validity of the argument based on chiral effective theory. 



As in [201, we have compiled data on the ratios A 30 (t)/Ai (t) and A 30 (t) / Ai (t) in Fig. 31 and Fig. 32 respectively, 
where we compare them to predictions based on experimental data for the nucleon form factors and PDFs in combi- 
nation with a model-dependent ansatz for the combined (x, independence of the GPDs, displayed by the error bands 
|93j . With increased statistics, even at the lightest pion mass there is good evidence that these ratios are generically 
not flat. A flat ratio would be a direct consequence of the factorization of the GPDs into an x-dependent function and 
a t-dependent function. In spite of a semi-quantitative agreement, we see that the predictions of Diehl et al. |93j are 
systematically steeper than the lattice data, in particular at small — t. This is most probably related to the fact the 
lattice data for the radii at the accessible pion masses are systematically below the experimental values. We further 
observe that the ratio of polarized GFFs is flatter than the unpolarized one (see for instance the m w ~ 600MeV graph, 
where the data is rather accurate). 

Figures 33 and 34 display respectively the ratios of GFFs C^o/^bo an d ^32/^301 both in the isovector and isosinglet 
channels. We have already remarked that C20 is consistent with zero in the isovector channel; however in the isosinglet 
channel, the ratio 4C2o/^20 is of the order of -1. This makes the £ dependence of H n=2 (^,t) and E n=2 (^,t) an order 
unity effect, see Eq. ([8]), which is of great interest due to its direct relation to the frequently discussed D-term [94]. For 
the next higher moment, n — 3, the data is rather noisy, but at = 495MeV and = 597MeV, it is nevertheless 
clear that 4A32/A30 is positive, and most likely of order +1. Thus the £ dependence of H n=3 (£,t) is also substantial, 
and goes in the other direction compared to the n = 2 sector. 



1. BChPT extrapolation of A^ d , B^ d , C^ d 

In this section we discuss the forward and small —t behavior of the generalized form factors A^q , -B20 arLC ^ ^20 
In particular A^q (t = 0) is the isovector momentum fraction (x) u —d- The main novelty as compared to |20) 
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FIG. 29: Generalized form factors of twist-two operators of dimension 3, 4 and 5, in the isovector and isosinglet channels. 
Notice in particular the flattening of the slope of the A n o GFF with increasing n. 
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that we now have sufficiently accurate data below m OT = 500MeV to test the applicability of covariant baryon chiral 
perturbation theory (BChPT) in that regime exclusively. In Ref. [20], the fit range extended up to = 700MeV. 
The 0(p 2 ) BChPT result [B] for the isovector GFFs {A,B,C}^ d {t) is 
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(61) 
(62) 

(63) 



Here fT d {m w ) 



h-A.B,c{ti m Tr) ancl j\ ( TO 7r) contain the non-analytic dependence on the pion mass and momentum 



transfer squared (see for instance Eq. (28), (40) and (41) of [71]), while A 2 q 



0,u-d 



and A 2 q d are low-energy constants. Similarly B, 



0,u-d 



= ^20 



we have included estimates of C(p 3 )-corrections in form of (5 t ^ u ~ d t), (5^ 7 " u ~ a m' 2 ), (5 l ( A u ~ a t) and (8™ 7 " M ~ a m' 2 ). The 
associated low-energy constants are treated as free parameters and may be obtained from a fit to the lattice data. 
Because of the small prefactor, the term cx hA[t,m n ) is of 0(1O~ 3 ) for < 700MeV, \t\ < lGeV 2 and therefore 
numerically negligible. Also, hB,c{t,m^) are very weakly dependent on t, so that the Ansatz for each form factor is 
practically linear in t. The parameter m^r in the denominator of Eqs. (62 1 and (63 1 is the nucleon mass in the chiral 
limit and was set to 890MeV, as given in Tab. |VIII| 

0.17 obtained from a heavy baryon chiral perturbation theory (HBChPT) fit to our 



'{t = 0, m w 



u Ji °20 = L/ 2 



A^ a (t = 0,m v = 0), A 



-d 



2(1 



(t = 0, m w = 0), and 



We use the value A 



0,u-d 
20 



TOtt < 500MeV lattice results for A 20 d (t = 0) = (x) / \u-Ad- This value is consistent with an earlier determination 



[95], which we used for a similar fit in [20]. The fit is displayed in Fig. 37 
Since the low energy constant A 



0,u-d 
20 



is a common parameter in the BChPT-formulae for the GFFs A 



t — d 
20 ' 



B. 



20 



and C 20 , we performed a simultaneous fit based on Eqs. (61 1 to (63 1 with a total of 9 (1 common and 



separate) 

fit parameters to over 83 lattice data points. The fit is displayed in the le ft panel of Fig. [35] and the results in the 

The x 2 /dof = 1.5 is good and gives us 
x) u -d — 0.1758(20), to be compared with 



XLIII 



chiral limit and at the physical pion mass are listed below in Tab. 
confidence that the fit works well. Among the most important results is 



i CTEQ6 



= 0.155(5). Our result is thus 10% to 15% higher than the phenomcnological value, a statistically significant 
difference. Compared to [20], it has increased from 0.157(10) by a little less than two standard deviations. The main 
reason is that A^ d (t = 0) at our lightest pion mass is bending down less than predicted by the chiral expansion 
compared to the heavier pion data points. 

We obtain B 2 ^ d (t = 0,mP hys ) = 0.293(12), a value compatible with our previous result 0.273(63) [20], but with a 



much reduced uncertainty. The quantity C 20 d (t = 0,mP hys ) 
and B 20 , but there is now a hint that it is in fact negative. 



-0.0157(86) still comes out much smaller than Ai 
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FIG. 31: Ratio of the A30 to the A\§ unpolarized GFFs, isovector and isosinglet. Disconnected contributions have been omitted. 
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FIG. 33: Ratios of n = 2 generalized form factors. 
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FIG. 34: Ratios of n = 3 generalized form factors. 
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FIG. 35: Simultaneous BChPT fit to isovector (left) and isosinglet (right) lattice data for < 580MeV. The dotted line is the 
heavy-baryon limit of the BChPT fit. A HBChPT fit to the lattice data for \t\ < 0.3GeV 2 and < 500MeV is shown by the 
dashed line. The phenomenological value of {x} u ±d (CTEQ6) is indicated by the cross. Disconnected contributions are omitted 
in the isosinglet case. Notice that B20 and C20 are displayed for t « -0.24GeV 2 . At m 2 w 0.12GeV 2 , the larger volume (28 3 ) 
is the filled symbol, the smaller volume (20 3 ) the open symbol; for all other pion masses, the volume is 20 3 . The x 2 /d.o.f. is 
1.5 and 2.8 respectively in the isovector and isosinglet cases. 



To study the difference between HBChPT and BChPT, we take the heavy baryon limit of BChPT while keeping 
the same values of the fit parameters, and obtained the dotted line in Fig. |35| This curve only overlaps with the 
BChPT curve for m, < m^ hys and drops off sharply for m w > mj tys , indicating the quantitative importance of the 
truncated terms when using the coefficients from the BChPT fit. The dashed curve in Fig. [35] shows the result of 
fitting our lattice data for \t\ < 0.3GeV 2 and < 500MeV directly with the HBChPT expression, and indicates 
that the latter describes the behavior of our lattice data over a significantly smaller range of pion masses than the 
BChPT expression. 

As already mentioned above, one of the main achievements of this work is the substantial reduction of the statistical 
uncertainties of the data points at = 356, 495 and 597MeV compared to our previous work |20j . and the inclusion 
of an additional ensemble at = 293MeV. This allows us in particular to lower the cut in in the chiral fits from 
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FIG. 37: The extrapolation of the polarized momentum fraction using HBChPT. The value obtained in the chiral limit is 
0.165(8), the experimental value is from HERMES [96]. 



700MeV to 600MeV. However, it is also important to note that a full BChPT analysis of ^207 -^20 an d C20 consistently 
including all terms of 0(p 3 ) is still not available. In view of this, it is interesting to perform a first check of the stability 
and potential uncertainty of the chiral extrapolations by repeating the fit for different maximal values of the included 
pion masses. Figure 38 shows a comparison of the BChPT extrapolations of from fits to the lattice data in the 



regions m n < 400, 500 and 600MeV. Most importantly, we find that the error bands from all three extrapolations 
are consistent and do overlap in the region below m\ ~ 0.25MeV 2 down to the chiral limit. Apparently, the bending 
towards the physical point is not overly sensitive to the large pion mass region, where 0(p 3 ) corrections would have 
the strongest impact. A quantitative estimate of the uncertainties of the chiral extrapolation must therefore be based 
on an improved ChPT analysis (e.g. including higher order effects) and higher precision lattice data at even lower 
pion masses, which is beyond the scope of the present work. 
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FIG. 38: A study of the stability of the BChPT extrapolation of at t = OGeV 2 based on Eqs. (|6Tj) , Q, and {63|. The 

three different error bands represent chiral fits to lattice results for pion masses in the regions < 400, 500 and 600MeV, 
respectively. 



2. BChPT extrapolation of A 



u-\-d r>u-\-d /~iu-\-d 
20 '-"20 i °20 



In this section we treat the GFFs in the isosinglet channel. The quantity A%q (t = 0) = {x) u +d is not only an 
important hadron structure observable on its own but in addition plays an important role for the computation of the 
total angular momentum contribution of quarks to the nucleon spin, J u + d = l/2(j4.2o~ d (0) + B^q (0)), a discussion 
that we postpone to Sec. IV E The combined (t, m^)-dependence in BChPT is given by [73] : 
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(66) 



where A% u+d = Aft 



(t = 0,m w = 0), and f^ +d (m- K ) and hA(t,m v ) contain the non-analytic dependence on the 



pion mass and momentum transfer squared. Also, B 



0,u+d rjM+d 



20 



20 



(t = Q,m n = 0), and the terms AB 



?t,u+d 



y0,U+d 



J B 

0,7717,- = 0), and the term AC 20 



20, 

-yu+d i 



t 



and 6™ w are of 0(p 3 ) and represent only a part of the full 0(p 3 ) contribution. Similarly, C^q 

d proportional to the pion momentum fraction carried by quarks in the chiral limit 



is a part of the full 0(p 3 )-corrections [74]. As in the isovector case, the parameter tojv in the denominator of 



Eqs. ([65]) and ([66]) was set to 890MeV. 

Thp constants A^,u+d At,u+d R 0,u+d <-t,u+d , c 
liic constants ^20 ' i , or alio o 



1 2(l 



'20 



B 



B 



may be obtained from a fit to the lattice data. In 



this counting scheme, contributions from operator insertions in the pion line, proportional to (x)^', d , are of order 



0(p 3 ). Counter terms of the form S t A u+d t and S^' u+d 777.2 g rst a pp ear a ^ 0(p A ). Since in [5D] the term AB^^it, m^) 
was found to lead to unstable fits, we perform a fit dropping this contribution but keeping the counter terms oc t 
and oc 7772. We have also included the formally higher order counter terms SQ U+d t and <5™"' u+<i 777^ in the fit to our 
lattice data. Furthermore, in order to see if such contributions could be relevant for the pion masses and values of 
the momentum transfer squared accessible in our calculation, we also include the estimate of the 0(p 3 )-contribution 
AA^q d provided in [71] in the fit to the lattice data points. 



Similar to the isovector case discussed in the previous sections, the low energy constant A. 



0,u+d 
20 



is a common 
— 0.5 as 



parameter in the chiral extrapolation formulae for the isosing let GFFs A£ a , B%+ a and C 2 " + • Using 
an input parameter J97H1Q0J , we performed a simultaneous fit to 83 lattice data points for these three GFFs, based 
on Eqs. (64), (65) and (66), with 1 common and 8 separate low energy constants as fit parameters. The results are 
also summarized in Tab. IXLlTTl 
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TABLE XLIII: Chirally extrapolated GFFs A 20 , B20 , C20 using BChPT. 
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FIG. 39: Pion mass and momentum transfer dependence of £20- At m n = 356MeV, the open symbol corresponds to 20 , the 
filled symbol to 28 3 . At all o ther pion masses, the volume is 20 3 . The extrapolation in m v at t = can be compared to the 



extrapolation of A^ d in Fig. 35 



The result x 2 /d.o.i. for this fit is 2.8. The result in the chiral limit is A, 



0,u+d 

20 



0.5567(45), and 



lu+d 



0, mP hys ) = 0.5534(43) at the physical point. Incidentally, this value is in very good agreement with phenomenological 
results from CTEQ and MRST [M] parameterizations, (x)^f T2001 = 0.538(22) and (x)^ Q6 = 0.537(22). A 



variation of the input parameter 



>7T,0 

lu+d 



by ±10% only leads to a small change in A 



o,»- 

20 



u+d 

(t = 0) of a few percent. 



The results of the fit are shown in Fig. 35 As already noted in [20], the slight upward bending in Fig. 35 at low 



TOtt, is due to the O(p a )-contribution AA^ d . The fact that the lightest pion mass data point does not exhibit this 
upward bending is partly responsible for the large value of the \ 2 - Further calculations are required to see whether, 
for example, this point is pushed downward by finite-size effects. The inclusion of contributions from disconnected 
diagrams could also lead to a different m n dependence. 



In Fig. 36 the simultaneous dependence of C%q on t and is displayed. The error bars of the lattice data 



points are illustrated by the stretched cuboids. The two surfaces represent the BChPT discussed in this section with 
its uncertainty. The figure illustrates that Cjq becomes relatively large near the origin, which implies that the £ 
dependence of ff" =2 (£,i) and E n=2 (£,t) is far from being negligible. 



3. HBChPT extrapolation of E^ d and M 2 " +d 

In heavy baryon chiral perturbation theory (HBChPT) [89, 90 to 0(p 2 ), the combined (i, m 7r )-dependence of the 
GFF-combination E%+ d (t) = Al+ d (t)+t/(Am 2 N )B%+ d {t) is quite different from that of M^ d (t) = A%+ d (t) + B%+ d (t) , 
which in the forward limit is equal to two times the total quark contribution to the nucleon spin, 2J q — (t = 0). 

The notation is chosen by analogy with the Sachs form factors in the n = 1 sector, E being the analog of Ge and 
M of Gm- Our results for B%+ d are in units of the quark-mass dependent nucleon mass; to be consistent with the 
conventions of the chiral expansion Eq. (68 1 below, we multiply the B^q^ data by m^ lys /r7ijv(ni 7r ) before performing 



the fit. In )20j . this conversion had been omitted. 

While at this order M^ +<i shows a non-analytic dependence on t and as discussed below, E"^ d is constant up 
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to analytic tree-level contributions, 
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(67) 



A fit to our lattice results based on Eq. ( |67| ) is shown in Fig. 39 The linear dependence of E 2 o on t and works 
well within the fit range, while the higher mass points do not quite lie on a common smooth curve. In contrast to 
the covariant approach, the functional form of E 20 (t — 0) = A 20 (t = 0) does not exhibit a term that could lead to an 
upward bending as seen on the right panel of Fig. |35| 
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FIG. 41: Pion and momentum transfer dependence of M20. At = 356MeV, the open symbol corresponds to 20 3 , the filled 
symbol to 28 3 . At all other pion masses, the volume is 20 3 . 



The pion mass dependence of M%Q~ d (t) for non-zero t is given by [8"9l |9"U] 



0,u+d 
20 



1 - 



39A m l i n f m 



M. 



;2,tt) 



Mm^.u+d 2 
20 "V 



M, 



t,u+d 



20 



(68) 



with new counter terms M^ t 7 " u+d and M^q^- The non-analytic dependence on t and in M^ 2 ' (t, m w ) results 
from pion-operator insertions and is directly proportional to the (isosinglet) momentum fraction of quarks in the pion 



in the chiral limit, {x)^ +d . We use (x)^ +d = 0.5 for the fit. The results of chiral fits based on Eq. (68) are presented 
in Fig. 41 We note that the data for M2o(t = — 0.24GeV 2 , m w ) has a tendency of bending downward as a function of 
Wtt, and the fit is able to describe this behavior. We find M^ +d (i = 0, mf ys ) = 0.528(11). We will use this result in 
Sec. EH 



4. HBChPT extrapolation of B%+ d 



Since the total anomalous gravitomagnetic moment of quarks and gluons in the nucleon has to vanish, Y^ q B2o(t — 
0) = 0, an interesting question is whether the individual quark and gluon contributions to B20 are separately zero or 
very small. The GFF B%Q d can be written as a linear combination of Eqs. (67) and (68). A separate fit to the data 



with fixed E. 



0,u+d 
20 



0.522 gives B%+ d (t = 0,mP hys ) = 0.015(11) which is compatible with the fits to M 2 " 4 



d and E^ d 



above that in combination give (M - E)^ d {t = 0,mP hys ) = 0.003(12) 



Although the absolute value of B£o~ d (i 



0) 

is again rather small, we note that the sign is different from that found in Sec. IV D 2 based on the BChPT fit, 
where it was negative. A more accurate calculation and the full 0(p 3 ) BChPT expression will help resolve the sign 
0). The simultaneous dependence of B^ 6, on and t is shown in Fig. 40 



rgu+d 
n 20 



We recall that the sum of 



and the corresponding quantity for gluons vanishes in the forward direction. The very non-trivial interplay of 
the t and dependence results in a small value of BV^ d near the origin of the (t, m^) plane, suggesting also a small 
value of -Bf (i = 0) in the real world. However, near the chiral limit the form factor -B^o becomes sizable, about -0.1 
at -t « 0.24GeV 2 . 



5. HBChPT extrapolation of C%o 
At order 0(p 2 ), the pion mass dependence of the GFF C^ d (t) is given by [88H90] 

C%+%m n ) = i _ f/ | 4m ^ ) |c 2 ° " +d + E£'*\t,m*) + E^\t,m,) + C m ^ d m\ + tj , (69) 

where C^o**^ = C^Q^it = 0,™^ = 0). The terms E% ' (t,m n ) and E% (t, m n ) contain non-analytic terms in t 
and to^ that come from insertions of pion operators proportional to (x)^ d . Additionally, E^"'^ (t,m n ) depends on 
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FIG. 42: Projections of C2 +d — as depicted in Fig. 



fit appears as a shaded band. 
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on a constant — t plane and on a constant m% plane. The HBChPT 



the low energy constants c%, C2 and C3, which are set to the same values as in |20j . To be consistent with the chiral 
expansion Eq. ( |69[ ), we rescale the C^ d lattice data by m^ lys /mjv(m 7r ) before performing the fit. We find a value 

(t = 0,mP hys ) = -0.325(14) at the physical pion mass. The 



of C, 



0,u+d 
20 



-0.407(14) in the chiral limit, and C 20 



u+d 



simultaneous dependence of CJo on m n and t is shown in Fig. 40 along with the HBChPT fit. Comparing the figure 
with Fig. [36j one sees that the fits differ significantly at larger values of m n and — t. The result is also shown as a 
function of the pion mass squared for fixed t « — 0.24GeV 2 in Fig. 
mass of 356MeV is presented in the right panel. 
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(left panel), and the i-dependence at a pion 



E. Quark contributions to the proton spin 



In the following, we will present and discuss our results for the quark spin and orbital angular momentum contri- 
butions (OAM) to the proton spin. We first remind the reader that the form factors of the energy momentum tensor, 
which in this section we denote by A q ' 9 (t) = A^it) and B q,9 (t) = i3| 9 (i), in the forward limit t — give direct 
access to the quark and gluon angular momenta in Ji's nucleon spin sum rule, 



q=u,d,... K q—u,d,. 



A q (Q) + B«(0) + A 9 (0) + B 9 {0) 



(70) 



Furthermore, the quark angular momentum may be decomposed in terms of the quark spin, AS 9 , and orbital angular 
momentum, L q , such that 



E 

q=u,d,. 



AS 9 + L q ]+ J 9 



(71) 



It is important to note that the decompositions in Eqs. |70j) and ( 71 ) are fully gauge-invariant, and that the individual 
terms will in general be renormalization scale and scheme dependent. Since the momentum fractions carried by quarks 
and gluons have to add up to one, i.e. the total nucleon momentum, 



1= A q (0)+A 9 (0)= J2 <^ + ^> 9 > 

q—u : d^... q—u,d,... 

one also finds that the anomalous gravitomagnetic moments, B q ' 9 (0), have to cancel exactly in sum, 

0= J2 B q (0) + B 9 (0). 

q=u,d,... 



(72) 



(73) 



It is a prominent goal of future lattice hadron structure calculations to study the above sum rules in great detail. 
Since gluonic observables suffer in general from very low signal-to-noise ratios and have so far not been studied on 
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the lattice with sufficient precision, for the moment we will have to concentrate on the (connected) contributions from 
up- and down qu arks. 



In Sec. 



IVE1 



below, we will begin with a discussion of J q based on our results for the GFFs A^q (t) and B^q (t). 



There, the results of the covariant and heavy baryon chiral extrapolations of Sees. IV D 1 IV D 2 and |IV D"3| will be 
supplemented and compared with an extrapolation of J u + d using a ChPT-formalism that includes explicitly the A- 
resonance as an additional degree of freedom. A decomposition of J q in quark spin and OAM contributions, together 
with corresponding chiral extrapolations, will be presented in Sec. |IV E 2| 



Quark angular momentum J 



From the covariant baryon chiral perturbation theory (BChPT) extrapolation in Sec. IV D 2 of the isosinglet GFFs 
A2o(t) and B2o(t), we find a value of JgchPT ~ 0.238(8) for the total quark angular momentum contribution at the 
physical pion mass. This corresponds to ~ 48% of the total nucleon spin S = 1/2, which is somewhat larger than our 
result in [50], although the difference is clearly not significant within statistical errors. We note again that these values 
have to be considered with some caution, since contributions from disconnected diagrams have not been included. 

A result that is more accurate regarding the systematics can be given for the iso vector, u — d, channel where 



IVD1 



we obtain J 



i—d 
BChPT 



0.234(6) 



disconnected diagrams cancel out exactly. From the chiral extrapolations in Sec. 

at mf ys . That this value is so close to the u + d-quark angular momentum already points to a small contribution from 
down quarks. Indeed, combining the isovector and isosinglet results, we find that the up quarks carry a substantial 
amount of angular momentum, JgchPT ~ 0.236(6), while the contribution from dow n quark s is very small and even 



IV E 2 below that the smallness 



negligible within the small statistical errors, J%ChPT = 0.0018(37). We will see in Sec. 
of J d can be traced back to a remarkably precise cancellati on b etween spin and orbital angular momentum of quarks. 

showing J u d 



The results we have just discussed are illustrated in Fig. 
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as a function of m 2 , together with the 



corresponding BChPT extrapolations indicated by the error bands. We note that the lattice data points were obtained 



from separate dipole extrapolations of the GFFs -B^W to t = 
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FIG. 43: Chiral extrapolations of J"' using BChPT. Note that the displayed lattice data points were not directly employed 
in the chiral fits. Details are given in the text. 

In order to study possible systematic uncertainties in the chiral extrapolation of quark orbital angular momentum 
just described, we now briefly discuss a completely different approach to the chiral extrapolation of J u+d . This will be 
based on HBChPT including the A resonance as an explicit degree of freedom, predicting at leading-one-loop order 
a pion mass dependence of J u+d of the form [102] 



J u+d (m ) 

HBChPT+A v'"^ 



I u+d (m ) 

HBChPTV "V 



iJVA 



'9(4tt/ t ) 



B 



0,u+d 
'20 



(ml -2A 2 )ln 



3(x 



7T,0 

u-\-d 



Y qA 



2A^A 2 - ml In 



+ VA 2 -ml J J 
and cjnNA is the pion-nucleon-A coupling. The m 



(74) 



where A = m& — m^ denotes the A-nucleon mass difference, 
dependent ^hbciiPT m ^1- (^4| corresponds to the HBChPT result without explicit A intermediate states as obtained 
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FIG. 44: Chiral extrapolation of J u+d using HBChPT including the A resonance, Eq. (74 1. The fit and error bands are 
explained in the text. 



from Eq. ( 66 1 for t = and is given by 



I u+d (m ) 

HBChPT V "V 



D 



0,u+d 

ha 



I u+d 



(A + B 



0,u+d 

ho 



\ 9A m l 

J i^hf 



111 



A 2 

X 



J 



(75) 



Since the GFF B 20 (t) cannot be extracted directly at t = 0, we have first performed separate dipole extrapolations 
of B% (t) and B$ (t) to t = 0, and combined this with our values for (x) u+d = A%+ d (0) to obtain J u+d = (A%+ d (0) + 
B 2Q ^ d (0))/2. The resulting lattice data points, including the full jackknife errors from the extrapolations of the 



B2o(t) to the forward limit, are displayed in Fig. 44 Chiral fits based on Eq. (74), with the three free parameters 
b qN = (A + B) 2 ^ +d , b qA and J™*> u + d , to the data with < 600MeV and < 700MeV are represented by the 
shaded error band and the curves (representing the upper and lower bounds of an error band) respe ctivel y. In both 
cases, we have fixed A = 0.3GeV and used the large- N c relation g^NA — 3/(2 3 ' 2 )gA as given in Tab. 



The fit to our lattice results with < 600 MeV gives (A + B 



ru-\-d 
HBChPT+A 



20 



VIII 



IqN 



0.514(41), b qA = 0.486(55) and 



(mP hys ) = 0.245(30) at the physical pion mass. Including the data point at = 687MeV in the fit, 



we find consistent values with somewhat smaller errors, (A + B)X + - b qN = 0.546(24) and b qA = 0.449(39) and 
^HBChPT+A( m £ hys ) ~ 0.226(22). It is encouraging to see that these values fully agree within statistical errors with the 
results from the global simultaneous BChPT extrapolations of the GFFs A 2 o(t), B 2 o(t) and C 2 o(t) discussed above. 



2. Quark spin and orbital angular momentum contributions 

For a consistent decomposition of the quark angular momentum, J 9 , into quark spin, AS 9 , and orbital angular 
momentum, L 9 , contributions, we need in addition lattice results for A^Q d (t = 0) = AT, u+d and A^Q d (t = 0) = 
AE"- rf . 



Our lattice data for AT, u+d /2 is displayed in Fig. 45 together with a 2-parameter HBChPT-fit represented by the 



upper shaded error band. The chiral extrapolation leads to a value of Af+ d /2(t = 0) = AS u+d /2 = 0.208(10) at the 
physical pion mass, perfectly matching the recent results from HERMES [96] indicated by the cross. However, since 
this is a leading 1-loop HBChPT fit at comparatively large pion masses, the agreement with the experimental value 
should be considered with great caution and seen as indicative. Combining this with the results from the previous 
section for J u + d and the corresponding BChPT extrapolation, we find a remarkably small quark orbital angular 
momentum L u+d = J u + d — AY, u+d /2 contribution to the nucleon spin for a wide range of pion masses, as indicated 
by the filled diamonds and the lower error band in Fig. |45[ From the combined covariant and heavy baryon chiral 
extrapolations, we obtain a value of L u+d — 0.030(12) at the physical pion mass. 

Superficially seen, such a small OAM contribution from u + d quarks of only w 6% to the nucleon spin is in 
clear conflict with general expectations from relativistic quark models, which suggest that L u+d = 30 — 40% of 1/2. 



3 In the notation of [102], b qN = {A + B)°^ +d = M°> u+d . 
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Moreover, from quite general arguments, e.g. based on light cone wave function representations of hadrons, substantial 
quark orbital motion is essential for the Pauli form factor Fi to be non- vanishing in general, and also for the formation 
of azimuthal single spin asymmetries in semi-inclusive deep inelastic scattering related to, e.g., the Sivers effect [1031 - 
1105] . As we will see in the following, these apparent inconsistencies may be explained by studying on the one hand the 
renormalization scale dependence of quark OAM, and on the other the contributions from individual quark flavors. 
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FIG. 45: Total quark spin and orbital angular momentum 
contributions to the spin of the proton. The cross represents 
the value from the HERMES 2007 measurement [92]. The 
error bands are explained in the text. Disconnected contri- 
butions are not included. 
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FIG. 46: Quark spin and orbital angular momentum con- 
tributions to the spin of the proton for up and down quarks. 
Filled and open squares denote AE u /2 and AE d /2, and filled 
and open diamonds denote L u and L d , respectively. The 
crosses represent the values from the HERMES 2007 mea- 
surement [96]. The error bands are explained in the text. 
Disconnected contributions are not included. 



We begin with the latter by noting that a study of the separate up- and down-quark OAM contributions requires in 
addition knowledge of the spin and angular momentum in the isovector, u — d, channel. To this end, we identify the 

AZ u ~ d , which we have 

cf. Fig. [9J Using our lattice data for g^, A£ u+d , and (A + -B)^**, we have computed 
the ind ividual spin and OAM contributions from up- and down quarks to the nucleon spin, which are displayed in 
Fig 



u — d quark spin contribution with the axial vector coupling constant, gA = A\ d {t = 0) 



discussed above in Sec. 



IV A 



4G 



as functions of . We stress that the chiral extrapolations represented by the shaded error bands in Fig. 



4G 



were not obtained from direct fits to the shown lattice data points, but are the result of combining the super-jackknifc 
error bands of the heavy baryon, SSE, and covariant baryon chiral extrapolations in Figs. [45] |9j and [35] The very 
good overlap of the bands with the data points may be seen as a first consistency-check of our approach, in particular 
with respect to the different types of extrapolations in the squared momentum transfer and the pion mass that we 
have employed throughout this work. 

Most remarkably, Fig. [46] shows that the individual up- and down-quark OAM contributions are sizable and of 
similar magnitude over a wide range of pion mas ses , but opposite in sign, and therefore only cancel to a large extent 
ju+d ag a i reac |y observed in Fig. 45 At the physical pion mass, we find from the chiral extrapolations 



in the sum, L u+ s» 0, as already observed in Fig. 

that \L U \ « \L d \ « 33% of 1/2. A more accurateresult can be given for the isovector channel where disconnected 
contributions cancel out, L u ~ d = -0.379(71). 

These observations may be seen in analogy to a corresponding analysis of the nucleon anomalous magnetic moment, 
k = F 2 (0) = £?io(0): Although the proton and neutron anomalous magnetic moments are sizable (and, as noted in 
[1031 1106] , related to non-zero quark orbital motion), k p = 1.79 and n n = —1.91, they largely cancel in the sum, 



3kP 



+n _ 



-0.36, see also Sec. 



IVB6 



However, as for L u+d , this doesn't imply that the orbital motion of 
individual quarks in the nucleon is necessarily small. 

Furthermore, from a direct comparison of Figs. [43] and [46] wc find that the smallness of the angular momentum of 
down-quarks can be seen as the result of another remarkable cancellation, in this case between spin and OAM. While 
both types of contributions for the down-quarks are similarly large in magnitude, |A£ d |/2 w \L d \ w 20 — 30% of 1/2, 
they are again of opposite sign and hence cancel out in the sum, J d — AE d /2 + L d w 0. 

Even though the small value for L u+d in our lattice calculation can be understood to arise from a cancellation 
between the different quark flavors, it would still be important to understand the striking discrepancy with results 
from relativistic quark model calculations, L u+d ps 30 — 40%, as noted above. In principle, this could be related to 
a severe deficiency of the model, or unexpectedly large, systematic uncertainties in the lattice calculation. A more 
likely solution has been proposed in [107[ 1108] . just by noting that quark model calculations generically correspond 
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to a low, hadronic scale [i <C lGeV, while the resolution scale chosen here is = 4GeV 2 . Since quark OAM is 
not conserved under (leading and higher order) evolution [109] , it is therefore in general pointless to directly compare 
the lattice results with the model expectations at the different scales. For a more sensible comparison, one might 
instead try to evolve the model results up to the lattice scale — or, alternatively, transform the lattice values down 
to the hadronic scale. In order to show that the different scales might be in principle responsible for the observed 
discrepancy, we follow Refs. [41, 108 and consider the LO scale dependence of L u ~ d , given by [109] 



L u - d (t) = 



f>0 V 



L u - d (t ) + -AS" 



1 



- -AS 1 
2 



(76) 



where i( ) — ^(^fo)/ ^QCd) with an initial scale denoted by [1q, and where the isovector spin contribution, AS" d /2 = 
gA/2, is exactly conserved under QCD evolution. Working in the isovector channel allows us to avoid issues related 
to disconnected diagrams (that were not included in our lattice calculation), and in particular the mixing with gluon 
operators under evolution. Figure 47 shows the result of the LO evolution according to Eq. ( 76 ) from an initial scale 
of no = 2GeV, with lattice starting values L u ~ d = —0.38 and AS"~ rf /2 = 0.61, down to very low (model) scales, as 
indicated by the left pointing arrow. 




/j 2 [GeV 2 ] 

FIG. 47: Leading order evolution of lattice results for L u ~ d (solid line) and AE"~ d (dashed line) from = 4GeV 2 down to 
low hadronic scales. 



Although the application of LO evolution equations at such low scales cannot be quantitatively trusted, we note at 
least that qualitatively a strong scale dependence sets in as /x — > Aqcd, eventually leading to a change of sign of L u ~ d . 
We note that similar observations have been made in the isosinglet channel |108) . even allowing for an approximate 
quantitative agreement of the lattice value for L u+d with the expectations from relativistic quark models. One should 
keep in mind, however, that the renormalization scheme (e.g. MS) of a model is in general indeterminate, which is 
a fundamental limitation in the comparison with results from the full theory. To sum up, the above exercise shows 
very clearly that a one-to-one confrontation of hadron structure observables from the lattice and in models must be 
considered with great caution. 



For convenience, we present in Tab. XLIV an overview of our results for the (connected) quark contributions to 
the proton spin, in the MS scheme at [i z — 4GeV 2 , at the physical pion mass, as obtained from the different chiral 
extrapolations discussed in this section and the previous sections. In the table, we also include estimates of systematic 



uncertainties due to the renormalization of the 1-derivative operators, as discussed in Sec. HID At this point we note 
again that no phenomenological values for AS = (1)a<3j ( x )q an d (%}Aq have been used in the chiral fits, and that we 
have so far only included contributions from connected diagrams in the lattice calculations. 

As a final note, we look forward to different angles of attack from the experimental and phenomenological sides 
becoming available. These will supplement the currently available DVCS measurements at JLab and HERMES and 
help to narrow down the statistical and systematic uncertainties [1101 1111] . Stronger constraints might come from 
experimental data on, e.g., the transverse target-spin asymmetry Ajjt for the electroproduction of vector mesons, 
which is in particular sensitive to the GPD E. In this respect, it is interesting to note that the (model dependent) 
phenomenological analysis in [112] , which is compatible with recent data from HERMES [113] and preliminary results 
from COMPASS (114] on Aut, has found a similar pattern of valence quark (orbital) angular momentum contributions 
as compared to our result in Tab. |XLIV| 
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V. SUMMARY AND CONCLUSIONS 



We have performed a comprehensive lattice QCD study of the observables characterizing the structure of the nucleon. 
The moments of forward parton distributions, the electromagnetic and axial form factors, and a set of generalized 
form factors, both for the isovector and the isoscalar flavor combinations were calculated. We have applied various 
kinds of chiral expansion schemes to the data and compared our results to experimental values. Determining a wide 
set of nucleon properties simultaneously allows us to assess the range of applicability of these schemes, and allows us 
to gauge how sensitive the nucleon properties are to changes in the up and down quark masses. 

The picture that emerges is that a nucleon composed of quarks corresponding to pion masses in the range 290MeV < 
to^ < 760MeV is significantly more compact than the nucleon in nature, when probed by a local vector current, as 
realized experimentally through photon-nucleon interaction. Both the Dirac and Pauli radii are less than a factor 
2/3 of their physical values, cf. Figs. 10 and 14 A similar observation holds for the axial radius, see Tab. |XXXIV 



consequently, a nucleon composed of heavier up and down quarks would also appear much smaller than the physical 
nucleon if probed by a W boson. When performing fits to our data using the small scale expansion (SSE) and 
covariant baryon chiral perturbation theory (BChPT) schemes we find that the pion mass dependence of our data is 
weaker than what we expect from those chiral expansions. Still, our findings are qualitatively compatible with the 
picture of a nucleon growing substantially in size when approaching the chiral limit. We also find that the Dirac 
form factor, at = 356MeV and for 0.3 < Q 2 /GeV 2 < 0.7, exhibits a statistically significant finite-size effect in the 
femto-universe in which we study the nucleon, see Fig. [H] The latter appears to be 'squeezed' when observed on the 
hypertorus of size (2.5fm) 3 . For this reason the chiral schemes might well be applicable to the pion mass regime which 
we probe if we had quantitative control of finite size effects. Hence, an important task for future lattice calculations 
will be to determine at what level such finite-size effects are influencing the extraction of the nucleon's infinite-volume 
properties. For the axial charge, we find the finite-size dependence to be at most (4 — 5)% at m T = 356MeV, cf. Fig. [3] 
but the effect is expected to grow in the chiral regime. 

We find the axial charge of the nucleon to be (8 — 10)% lower than in the physical world, and its pion mass 
dependence is very weak in the explored range, see Fig. [9] The latter observation is in qualitative agreement with 
the prediction of chiral effective theory, and continues up to pion masses far outside its range of validity. It is an 
important challenge to show that gA really rises by (8 — 10)% as is decreased from 300MeV down to its physical 
value. 

We have also presented extensive data on the generalized form factors (GFFs), which are related to the s-moments 
of the generalized parton distributions (GPDs). The higher the spin of the twist-two operator, the smaller the radii of 
the associated GFFs. This observation holds throughout the explored pion mass range, see Fig. [29] and is in accord 
with the expectation of a decreasing transverse size of the nucleon, as measured by the impact parameter of an active 
quark relative to the nucleon center of momentum, with increasing longitudinal momentum fraction x. Since only 
the vector and axial-vector radii (i.e. for the lowest moment n = 1) exhibit a statistically significant growth when the 
pion mass is reduced, cf. Fig. |30[ the hierarchy between these radii and those associated with higher spin operators 
increases towards the physical point. We have also found that the longitudinal momentum transfer £ dep end enc e of 
the GPDs is significant, as the ratios C^ d /A^+ d and A%+ d /A%+ d show quite clearly, cf. Eq. (|| and Figs. [33] and [34 
This is in fact consistent with large-7V c counting rules, which predict in particular that C^Q d ~ N c ~ A^ • 

The GFFs in the forward limit, namely the moments of parton distributions, exhibit a rather mild pion mass de- 
pendence in the accessible range of pion masses. Specifically, both the unpolarized and polarized isovector momentum 
fractions, shown in Figs. [35] and [37] respectively, will have to bend down as a function of beyond our range, if the 
lattice data is to make contact with experiment. 

Finally, we applied our results for the moments of the GPDs to the decomposition of the proton spin — the final 



results are listed in Tab. XLIV Thus, we extend our previous calculation in Ref. [20] down to = 290MeV, and 
increase its accuracy. The decomposition of quark angular momentum in terms of quark spin and quark orbital angular 
momentum is displayed in Fig. [46] The spin contributions are in agreement with the HERMES 2007 data while 
orbital contributions are our predictions. We remind the reader, however, that our calculation did not include the 
disconnected graphs which contribute to all but isovector quantities, and that the mixing of the isosinglet quark 
operators with the gluonic operators was not taken into account. These systematic uncertainties should be kept in 
mind in the interpretation of the non-isovector results. The isovector quantities however, such as ^A'E u ~ d and L u ~ d , 
do not suffer from these systematic uncertainties. One of the surprises, from the relativistic quark-model point of 
view, is the negative sign of L u ~ d . Various interpretations of this result have been proposed |108[ 1115] , As pointed 
out in Ref. [108] . L u ~ d evolves rapidly at low renormalization scales, because J u ~ d = L u ~ d + ^AY, u ~ d renormalizes 
multiplicatively and \AY, u ~ d is scale-invariant. In any case, lattice studies of the Ji sum rule, Eq. (70), have led to 
a renewed interest in the problem of 'decomposing the total angular momentum of an interacting multiconstituent 
system into contributions from various constituents' (cited from [115] : see also [116 ). 

Given that lattice QCD provides a systematically improvable way of solving QCD, there should be eventual agree- 
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ment between nucleon properties calculated in this framework and those measured by experiment. Thus, an extension 
of the present calculation down to = 200MeV is almost certain to uncover the dramatic effects predicted by chiral 
effective theory near the chiral limit. To be specific, the isovector Dirac radius (rf) and the Pauli radius (rf) can 
provide important benchmarks of our understanding of nucleon structure — they both need to grow particularly fast 
as m-n- decreases below 250MeV, if our lattice data are to make contact with the phenomenological values, see Figs. [12] 
and |16[ The calculation of these and similar observables would thus test our understanding of both chiral effective 
field theory and of lattice calculations as applied to nucleon structure. 

To achieve this goal we must be sure that statistical uncertainties, cut-off effects, and finite-size effects are well 
understood and under control. We believe that our current work provides an important step in this direction, as 
it demonstrates the improvement we have achieved in reducing statistical uncertainties and the application of chiral 
perturbation theory. In view of the success of our methods and techniques and given the recent increase in the 
computing power available to lattice QCD practitioners, we are optimistic that this program can be carried out 
successfully in the course of the next few years. 
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